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ABSTRACT 


M-Layer,  the  tropospheric  propagation  effect  prediction  program  by  NRaD  (formerly 
NOSC),  is  revised  for  greater  accuracy,  speed  and  stability.  This  is  achieved  through 
converting  the  extended  complex  number  representation  into  the  representation  by  the 
complex  exponent,  improving  the  accuracy  in  Airy  function  computation,  introducing  a 
new  mode  locating  algorithm  and  implementing  a  consistency  checking  procedure  for 
determining  the  proper  method  to  evaluate  the  height  gain  function.  The  revision  has 
been  documented  and  the  new  program  source  code  has  been  delivered  to  NRaD.  It  is 
recommended  that  the  mode  search  protocol,  not  just  the  mode  locating  algorithm 
introduced  in  this  revision,  be  completely  revised.  Unlike  the  current  approach  of 
blanketing  the  whole  possible  region  until  exhaustion,  modes  should  be  searched 
according  to  their  range  attenuation  rates  one  by  one  along  a  well  defined  path.  This 
should  result  in  a  faster  and  even  more  stable  program.  The  program  size  can  also  be 
reduced. 
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I.  INTRODUCTION 


M-Layer  is  a  FORTRAN  program  for  computing  the  propagation  factor  of  an 
electromagnetic  (EM)  wave  in  a  stratified  atmosphere.  It  is  desirable  to  extend  the 
capability  of  this  program  to  include  a  layer  of  random  medium  representing  the  air- 
ocean  interface  where  the  thickness  of  this  layer  cannot  be  ignored,  where  the  EM 
propagation  and  scattering  are  so  strongly  coupled  that  clutter  and  propagation 
effects  within  this  layer  cannot  be  dealt  with  separately,  and  where  the  grazing  angle 
of  the  EM  wave  incident  into  this  layer  is  so  small  that  the  curvature  of  the  earth 
cannot  be  neglected.  To  achieve  this  goal,  there  are  many  basic  theoretical  problems 
which  have  to  be  answered.  First  of  all,  the  effect  of  the  earth  curvature  in  this 
program  is  taken  care  of  through  the  classical  earth-tlattening  approximation  [Ref. 
l],but  the  result  [Ref.  2]  does  not  agree  with  the  more  recent  diffraction  theory  of 
Fock  [Ref.  3]  near  the  surface  of  the  earth.  Then  there  is  the  question  about  the 
better  method  to  model  the  atmospheric  refractive  index  profile,  either  piecewise 
linear  or  quadratic,  to  be  resolved  by  a  new  earth-flattening  approximation  under 
development  at  NFS.  The  new  approximation  w'ill  also  determine  the  functions  to  be 
used  for  the  representation  of  the  EM  fields  in  each  layer  through  uniform 
asymptotic  theories.  Within  some  proper  region,  these  new  functions  are  expected  to 
reduce  to  the  Airy  functions  utilized  by  M-Layer.  The  evolutionary  nature  of  this 
effort  prompted  this  review  to  improve  the  inner  workings  of  the  M-Layer  program. 
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In  particular,  the  subroutines  to  search  for  the  modes  and  those  for  evaluating  the 
Airy  functions  will  remain  as  an  important  part  of  a  program  investigating  questions 
about  EM  wave  propagation  by  solving  the  related  boundary  value  problem. 

It  can  never  be  overemphasized  that  a  boundary  value  problem  which  includes 
a  layer  of  random  medium  or  some  range  dependent  inhomogeneity,  set  up  according 
to  the  Maxwell  equations,  wilt  include  backscattering  in  its  solution.  This  is  in  sharp 
contrast  to  those  numerical  procedures  based  on  the  parabolic  approximation  to  the 
wave  equation  for  which  the  backscattering  is  completely  ignored. 

In  what  follows,  the  M-Layer  program  and  the  reasons  for  replacing  the 
extended  complex  numbers  with  their  complex  exponent  representations  are 
discussed,  together  with  some  other  problems  encountered  and  resolved  during  this 
investigation. 

A.  M-LAYER 

In  M-Layer,  the  index  of  refraction  of  the  atmosphere  is  assumed  to  be  height 
dependent  and  is  approximated  with  a  continuous  piecewise  linear  profile.  The 
classical  earth-flattening  approximation  is  utilized  to  allow  the  use  of  the  cylindrical 
coordinate  system  while  retaining  the  effect  of  the  curv-ature  of  the  earth.  This  is 
done  simply  by  substituting  the  index  of  refraction  with  the  modified  index  of 
refraction,  which  also  has  a  piecewise  linear  profile  [Ref  1]. 

The  source  of  the  EM  radiation  is  assumed  to  be  either  a  vertical  electric 
dipole  or  a  vertical  "magnetic  dipole’,  with  the  latter  providing  an  approximation  to 


the  radiation  of  a  horizontal  electric  dipole.  The  dipole  is  located  along  the  positive 
z-axis  of  the  cylindrical  coordinate  system  while  the  origin  is  sitting  on  the  ground. 
The  x-y  plane  is  the  "flattened"  earth  surface.  After  carrying  out  the  Hankel 
transform  along  the  radial  direction,  the  resulting  spectrum  of  the  Hertzian  dipole 
field  within  each  layer  of  a  linear  segment  of  the  modified  refractive  index  profile  is 
reduced  to  a  linear  combination  of  the  Airy  functions.  Specifically,  the  layers  are 
numbered  to  increase  with  height,  with  the  first  layer  being  the  one  right  above  the 
ground.  The  spectrum  of  the  Hertzian  dipole  field  is  proportional  to  the  product  of 
the  values,  at  the  transmitter  height  and  at  the  receiver  height  respectively,  of  the 
height-gain  function.  At  a  height  within  the  i-th  layer,  the  height-gain  function  is 
given  by  [Ref.  4]: 

/.(p,z)=6,(pMA,(p)fc,(g,)^fc,(9,)]  ,  (1) 

where  p  is  the  radial  component  of  the  propagation  vector  and  is  also  the  spectral 
variable  of  the  Hankel  transform;  hence  it  is  the  same  throughout  all  layers.  It  is  a 
complex  variable  whose  imaginary  part  represents  the  radial  attenuation  rate  of  the 
spectral  component  of  the  Hertzian  dipole  field.  Under  the  classical  earth-flattening 
approximation,  the  spectrum  of  the  Hertzian  dipole  field  contains  a  discrete  portion 
and  a  branch  cut.  The  discrete  spectrum  gives  rise  to  the  creeping  wave  modes 
diffracted  by  the  earth  surface  and  the  dielectric  waveguide  modes  supported  by  the 
layered  atmosphere.  The  contribution  from  the  branch  cut  is  usually  negligible, 
especially  for  the  field  in  the  shadow  of  the  earth.  The  M-Layer  program  locates  the 


discrete  spectrum  for  modes  having  a  radial  attenuation  rate  below  a  predetermined 
value.  Contributions  from  these  modes  determine  the  propagation  factor  of  the  wave. 

The  variable  ry,  in  the  i-th  layer  is  a  dimensionless  linear  function  of  height  z 
with  the  free  space  wavenumber  the  modified  index  of  refraction  w-  at  the  lower 
boundary  z  =  z,,  the  slope  of  the  modified  index  of  refraction  a, 72  and  p  as 
parameters: 


'1 


/ «.  \y 


2  /  N  P 

rtt:  +a.(z-z)-— 
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The  height  dependence  of  the  field  is  given  in  terms  of  the  functions  k^{q^)  and 
kiiOj),  w'hich  are  proportional  to  the  Airy  functions  and  Ai{—q,) 

respectively.  Of  these  two  functions,  at  a  height  so  large  that  ry,  is  large  and  positive, 
^i(ry,)  represents  a  downward  going  wave  and  e''^’^^^A'i(fy,)  +  A’2(ry, )  represents  an 
upward  going  wave.  The  coefficients  /t,  and  B-  are  determined  by  the  conditions  on 
the  continuity  of  the  Hertzian  dipole  field  and  its  derivative  across  layer  boundaries 
and  by  the  normalization  condition  that  the  integral  ,  :  the  square  of  the  height-gain 
function  over  all  height  equals  unity. 

To  fulfill  the  radiation  condition,  the  highest  layer  is  given  the  same  refractive 
index  as  the  free  space  above  it  and  only  the  outgoing  wave  is  allowed  within  this 
layer.  Below  the  "flattened"  earth  surface,  the  field  is  assumed  to  be  a  plane  wave 
propagating  downward.  Hence  only  the  normalization  factors  are  required  in  the 
highest  layer  and  in  the  ground.  By  assigning  B,  to  unity  in  the  highest  layer,  all  the 
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coefficients  /I,  and  can  be  determined,  according  to  the  boundary  conditions,  to 
within  a  multiplicative  factor  for  This  multiplicative  factor  is  then  deduced  from 
the  normalization  condition.  This  procedure  can  also  be  carried  out  from  the  ground 
level  up.  That  these  coefficients  can  be  computed  either  from  the  highest  level  down 
or  from  the  lowest  level  up  is  a  result  of  the  fact  that  p  belongs  to  the  discrete 
spectrum  of  the  Hertzian  dipole  field.  Consequently,  agreement  between  these  two 
ways  of  evaluating  the  and  coefficients  confirms  that  a  mode  has  been  located 
accurately. 

B.  EXTENDED  COMPLEX  NUMBER  REPRESENTATION 

The  discrete  spectrum  of  the  Hertzian  dipole  field  corresponds  to  the  zeroes 
of  the  modal  function  which  is  a  determinant  whose  elements  consist  of  k^{q^)  and 
AsCiy,)  at  the  layer  boundaries.  Numerically,  the  magnitude  of  this  modal  function 
causes  overflow  and  underflow  problems  as  k^Ul,)  or  Astry,)  becomes  exponentially 
large  or  small  for  complex  q-  values.  In  the  M-Layer  program,  to  overcome  this 
problem,  a  complex  number  is  written  as  a  scaled  number,  which  is  complex, 
multiplied  by  a  scaling  factor  which  is  an  integer  power  of  e,  the  base  of  natural 
logarithm.  This  integer  is  chosen  so  that  the  greater  of  the  absolute  values  of  the  real 
part  and  the  imaginary  part  of  the  scaled  number  lies  within  e-’.  A  complex  number 
written  in  this  form  is  called  an  extended  complex  number.  Multiplication  of  two 
extended  complex  numbers  requires  summing  the  two  integer  exponents  in  addition 
to  carrying  out  the  regular  complex  multiplication  of  the  scaled  numbers.  Addition 


of  two  such  numbers  is  achieved  through  the  use  of  an  addition  subroutine:  the  larger 
scaling  factor  is  factored  out  of  both  addends  before  they  are  combined.  The  scaling 
factor  is  adjusted  after  each  addition  and  after  a  sequence  of  multiplications  to  make 
sure  that  the  resulting  scaled  number  is  still  within  the  desired  range.  Addition  is 
troublesome  when  the  two  numbers  to  be  added  nearly  cancel  each  other.  Under  this 
circumstance,  the  scaling  factors  of  the  two  numbers  are  identical  and  both  the  real 
parts  and  the  imaginary  parts  of  the  scaled  numbers  are  almost  equal  with  opposite 
signs.  It  is  clear  that  the  real  part  and  the  imaginary  part  of  the  sum  lose  their 
accuracies  to  different  degrees;  hence  the  phase  angle  may  incur  substantial  error. 
To  remedy  this  situation,  interpolation  procedures  have  to  be  devised. 

As  two  complex  numbers  come  close  to  cancel  each  other,  they  must  be  out  of 
phase  by  almost  180  degrees.  By  factoring  out  the  square  root  of  their  product 
instead  of  the  scale  factor,  the  resulting  addends  become  reciprocal  to  each  other, 
both  lying  within  an  identical  small  angle  to,  and  on  the  same  side  of,  the  imaginary 
axis.  They  are  close  to  the  unit  circle,  but  one  is  on  the  inside  and  the  other  is  on  the 
outside.  Taking  out  further  a  phase  factor  of  -kH  after  writing  the  addends  in  their 
exponential  forms,  the  exponents  become  small  numbers  for  w'hich  Taylor  series 
expansion  of  the  exponential  function  converges  rapidly  and  can  be  used  for 
interpolating  the  sum  to  achieve  higher  accuracy.  Note  that  after  the  extra  phase 
factor  of  ■jr/2  is  removed  from  the  addends,  it  is  actually  the  difference  of  the 
resulting  two  reciprocals  which  is  computed.  This  procedure  effectively  picks  the 
direction  on  the  complex  plane  along  which  the  addends  are  almost  opposing  each 
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other  to  carry  out  their  cancellation.  The  resulting  sum  has  a  phase  angle  nearly 
perpendicular  to  this  chosen  direction. 

It  is  evident  that  the  representation  of  a  complex  number  by  its  complex 
exponent  of  base  e  provides  better  phase  accuracy  for  addition.  A  one-to-one 
correspondence  can  be  achieved  by  restricting  the  imaginary  part  of  this  exponent  to 
within  -X  and  x.  This  will  be  called  the  exponential  representation  or  the  complex 
exponent  representation  henceforth.  It  is  convenient  for  multiplication:  adding  the 
complex  exponents  of  the  two  factors  will  suffice.  Conversion  of  the  M-Layer 
program  from  the  extended  complex  number  to  the  complex  exponent  representation 
has  been  carried  out. 

C.  CONSISTENCY  CHECKING 

As  better  precision  is  achieved,  problems  with  the  mode  search  procedure  and 
the  evaluation  of  the  A^  and  S,  coefficients  become  severe.  They  are  thoroughly 
investigated  and  resolved.  For  mode  search,  although  the  division  of  the  region  of 
interest  into  "contour  rectangles"  and  further  into  square  "meshes"  and  the  search 
pattern  to  move  around  the  sides  of  a  "contour  rectangle"  to  find  and  follow  "phase 
lines"  into  it  are  kept,  the  basic  assumption  of  Shellman  and  Morfitt  [Ref.  5]  that 
both  the  real  and  the  imaginary  parts  of  the  modal  function  are  linear  along  every 
edge  of  a  mesh  square  is  completely  abandoned.  For  the  evaluation  of  the  /I,  and  fi, 
coefficients,  the  "test  for  evanescence"  conditions  have  been  removed.  A  consistency 
condition  to  determine  whether  to  evaluate  the  coefficients  from  the  ground  level  up 
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or  from  the  top  level  down  has  been  fomulated  and  incorporated  into  the  program. 
This  accomplishment  leads  to  the  relaxation  of  mode  locating  accuracy  requirement 
which,  combined  with  the  improved  precision  of  the  revised  program,  makes  the  first 
order  Newton-Raphson  iteration  unnecessary.  The  specific  changes  in  the  program 
and  the  resulting  gains  in  speed,  accuracy  and  execution  stability  are  discussed  in  the 
following  chapters.  Recommendation  to  completely  revise  the  mode  search  protocol 
to  do  without  the  "contour  rectangles"  is  also  provided. 
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II.  PROGRAM  REVISIONS 


M-Layer  is  structured  into  three  parts:  setup,  mode  search  and  propagation 
factor  evaluation.  The  main  input  is  the  modified  refractive  index  values  at  specified 
heights  so  that  a  piecewise  linear  profile  can  be  constructed.  If  the  mode  locations 
for  the  particular  profile  are  available  from  a  previous  run  of  the  program,  they  can 
also  be  included  in  the  input  and  the  mode  search  procedures  will  be  bypassed.  The 
various  ranges  and  transmitter  and  receiver  heights  for  which  propagation  factors  are 
desired  are  also  specified.  The  subroutine  WVGSTDIN  is  called  to  input  the 
information  from  an  ASCII  data  file.  The  program  then  computes  the  constants  to 
be  used  for  mode  search  and  propagation  factor  evaluation.  The  mode  search  is 
performed  with  the  subroutine  FNDMOD.  The  MODSUM  subroutine  is  then 
invoked  to  first  compute  the  and  S,  coefficients  as  explained  in  the  Introduction, 
then  compute  the  propagation  factor  and  the  propagation  loss.  The  complete 
program  structure  is  given  in  Figures  1  and  2.  There  are  several  other  subroutines 
which  are  not  included  in  these  and  other  figures,  such  as  DHORIZ  for  computing 
the  horizon  distance  between  a  transmitter  and  a  receiver  for  reference  purpose; 
CHKMOD,  a  maintenance  routine  for  removing  zero  from  reported  mode  locations 
by  older  versions  of  the  program;  or  A02H20,  a  routine  to  compute  the  atmospheric 
absorption  coefficient  due  to  oxygen  and  water  vapor.  They  will  not  be  discussed  as 
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they  do  not  contribute  directly  to  the  main  purpose  of  this  program  of  locating  the 
modes  and  computing  the  propagation  factor. 

The  program  structure  has  been  altered  as  shown  in  Figures  3  and  4.  Since  the 
/I,  and  B-  coefficients  have  to  be  evaluated  only  once,  they  are  now  obtained  through 
a  call  to  the  subroutine  AECOEF  directly  from  the  main  program  right  after  the 
modes  are  located.  Several  subroutines  are  dropped  in  this  revision  for  various 
reasons:  The  subroutines  NORME  and  NORMRE  are  eliminated  because  they  are 
no  longer  needed  due  to  the  change  in  complex  number  representation;  The 
subroutines  NOMSHX,  EDFDTX  and  DXDETR  are  not  used  because  the  modes 
are  now  located  with  adequate  precision  without  further  iteration;  The  subroutine 
ADDX  is  not  listed  separately  because  it  is  called  only  once  and  has  been  reduced 
to  only  a  few  lines  which  are  placed  where  the  subroutine  is  called  in  the  original 
program.  On  the  other  hand,  changes  in  the  mode  search  algorithm  require  the 
addition  of  two  new  subroutines:  SURFO  is  a  modified  and  simpler  version  of  SURF; 
ROOTS  replaces  QUAD.  Due  to  the  change  in  complex  number  representation,  all 
subroutines  listed  below  FNDMOD  and  MODSUM  have  been  revised,  including 
their  input/output  lists.  But  except  for  SURFO  and  ROOTS,  the  utilities  of  these 
subroutines  are  the  same  as  those  of  the  original  ones.  Descriptions  of  these 
subroutines  can  be  found  in  the  report  by  Yeoh  [Ref.  4]. 

The  most  significant  changes  have  been  made  in  XCADD,  XCDAIT  and 
XCDAIG  for  adopting  the  complex  exponent  representation  and  improving 
computation  speed  and  accuracy;  in  FZEROX,  FINDFX,  ROOTS  and  SURFO  for 


stabilizing  and  simplifying  the  mode  search  algorithm;  and  in  ABCOEF  for 
implementing  the  criteria  to  determine  the  reliable  manner  for  evaluating  the  A  -  and 
fiy  coefficients.  These  changes  are  discussed  in  the  sections  below.  The  source  code 
listings  of  the  completely  new  subroutines  XCADD  and  ROOTS  and  the  significantly 
revised  subroutines  FZEROX  and  ABCOEF,  which  are  compiled  with  Microsoft 
FORTRAN  version  5. 00, are  attached  as  Appendices  A  through  D.  Validation  of  the 
revised  program  has  been  carried  out  at  9.6  GHz  for  all  the  21  profiles  listed  in 
Yeoh  [Ref.  4]. 

A.  ADDITION  SUBROUTINE 

XCADD  is  the  subroutine  implementing  the  addition  of  complex  numbers 
under  the  representation  by  their  exponents.  Given  the  double  precision  complex 
numbers  Zj  and  Zi  as  the  exponents  of  the  addends,  this  subroutine  returns  the 
exponent  of  the  sum.  Since  a  double  precision  number  has  an  accuracy  of  53  bits,  if 
the  real  parts  of  z,  and  Z2  differ  by  more  than  53  bits,  the  exponent  of  their  sum  will 
simply  be  the  one  of  the  greater  real  part.  When  cancellation  becomes  serious,  the 
square  root  of  the  addends  is  factored  out  first.  Then  the  four-term  Taylor  series 
expansions  of  the  resulting  reciprocals  are  summed  up.  Since  the  leading  term  of  the 
sum  of  the  Taylor  series  is  a  good  estimate  of  the  sum  of  the  reciprocals  and  the 
relative  error  of  the  four-term  Taylor  series  sum  is  proportional  to  the  fourth  order 
of  this  leading  term,  the  threshold  for  invoking  this  interpolation  procedure  is  set  at 
the  highest  possible  value  of  2“*'^  allowed  under  double  precision.  Experimenting 
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with  this  procedure  shows  that  this  interpolation  improves  accuracy  as  long  as  the 
threshold  is  set  at  a  number  between  2“^"*  and 

B.  AIRY  FUNCTION  EVALUATION 

Similar  to  the  original  program,  the  evaluation  of  the  Airy  function  adopted  the 
algorithm  prescribed  by  Schulten,  et.  al.  [Ref.  6].  In  the  new  program,  changes  are 
made  to  follow  the  advice  of  Schulten,  et.  al.  concerning  the  region  within  which 
Taylor  series  expansion,  instead  of  the  faster  Gaussian  quadrature,  has  to  be  used  to 
achieve  double  precision  accuracy.  Other  changes  in  implementing  the  algorithm  are 
described  below. 

1.  XCDAIT 

Due  to  the  similarity  in  their  Taylor  series  coefficients,  the  Airy  function 
and  its  derivative  are  evaluated  within  a  single  loop.  The  relative  accuracy  of  the 
derivative  of  the  Airy  function  is  set  at  the  double  precision  limit  of  2”^“*. 

2.  XCDAIG 

Six  term  Gaussion  quadrature  is  used  for  evaluating  the  Airy  function  and 
its  derivative  outside  the  circle  of  radius  4.97 centered  at  (0.90, 2.80)  on  the  complex 
plane.  The  use  of  four-term  quadrature  outside  a  radius  of  15  from  the  origin 
suggested  by  Schulten,  et.  al.  is  not  adopted.  The  six-term  quadrature  in  this  range 
retains  a  higher  accuracy  while  overall  speed  improvement  by  using  both  the  four- 
term  and  the  six-term  quadrature  appears  to  be  minimal. 
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C.  MODE  LOCATING 


As  explained  in  the  Introduction,  the  modes  are  located  at  the  zeroes  of  the 
modal  function.  These  zeroes  are  located  on  the  upper  complex  qjj  plane.  Here  qjj 
is  the  value  of  qj  on  the  earth  surface,  which,  according  to  Eq.(2)  of  Chapter  1,  is  a 
linear  function  of  p“.  For  a  horizontally  propagating  mode,  p/k  is  close  to  unity.  The 
maximum  range  attenuation  rate  specified  for  the  desired  modes,  which  corresponds 
to  a  limit  on  the  imaginary  part  of  p,  determines  approximately  the  upper  bound  for 
the  imaginary  part  of  the  q^  complex  plane  to  be  searched  for  modes.  The  Shellman 
and  Morffit  mode  search  procedure  first  divides  the  search  region  horizontally  into 
"contour  rectangles"  each  of  which  spans  160  meshes  along  the  real  qn  direction.  A 
mesh  is  a  square  whose  size  is  an  adjustable  parameter  of  the  order  10""^  at  9.6GHz 
for  most  of  the  cases  considered  herein.  This  parameter  is  determined  by  the 
frequency  and  the  slope  of  the  modified  index  of  reflection  in  the  lowest  layer  of  the 
profile.  The  search  commences  at  the  top  left  corner  of  the  "contour  rectangle"  whose 
left  edge  has  a  real  coordinate  value  close  to  the  difference  of  the  real  parts  of  the 
qii  values  with  the  minimum  modified  index  of  refraction  and  the  index  near  the 
surface  substituted  into  Eq.(2)  of  Chapter  1.  After  the  search  over  the  initial 
rectangle  is  completed,  the  program  moves  to  search  the  next  rectangle  until  a 
specified  maximum  number  of  modes  are  found  or  a  specified  number  of  "contour 
rectangles"  have  been  searched. 

The  search  for  zeroes  makes  use  of  the  fact  that  a  real  function  changes  sign 
w'hen  it  crosses  a  simple  zero.  Since  a  zero  of  a  complex  valued  function  F(q)  is 
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where  both  its  real  part  and  imaginary  part  vanish,  a  necessary  condition  for  a  point 
to  be  a  zero  is  that  it  is  on  the  intersection  of  two  curves  defined  by  Im{F(q)}=0 
and  Re{F(q)}=0.  The  program  searches  around  a  "contour  rectangle"  for  a  sign 
change  in  Im{F(q)}  across  an  edge  of  a  mesh  bordering  the  side  of  the  "contour 
rectangle"  to  determine  that  a  line  of  Im{F(q)}=0  has  been  encountered.  The  search 
then  follows  this  line  into  the  meshes  within  the  "contour  rectangle',  checking  each 
mesh  to  see  if  a  cur\'e  Re{F(q))=0  enters  the  mesh  under  investigation.  All  these 
steps  make  use  only  of  the  assumption  that  the  zeroes  of  the  modal  function  are 
simple.  Once  both  the  cur\e  Im{F(q)}=0  and  the  curve  Re{F(q)}=0  are  determined 
to  be  present  within  a  mesh,  the  location  of  their  possible  interception  is  estimated. 
An  algorithm  for  this  estimate  is  required. 

Shellman  and  Morffit  [Ref.  5]  introduced  a  further  assumption  that  the 
functions  Re{F(q)}  and  Im{F(q)}  are  both  linear  along  the  edges  of  a  mesh.  Based 
on  this  assumption,  they  try  to  estimate  the  locations  where  the  curve  Im{F(q)}=0 
enters  and  leaves  a  mesh  square  and  the  location  of  q,^  if  a  curve  Re{F(q)}=0  also 
enters  the  same  mesh.  It  is  obvious  that  information  about  the  locations  w'here  the 
curves  enter  and  leave  the  mesh  square  is  not  essential.  Furthermore,  in  the  18  m 
duct  height  case,  the  scheme  causes  the  search  path  to  loop  around  four  contiguous 
meshes  until  the  search  is  broken  up  by  the  limit  on  the  number  of  meshes  to  be 
investigated.  Replacing  their  technique  requires  major  changes  in  the  subroutines 
involved.  A  new  subroutine  ROOTS  is  provided  to  estimate  the  location  of  the 
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intersection  of  the  curv'es  lm{F(q)}=0  and  Re{F(q)}=0.  These  changes  eliminate 
the  looping  problem. 

Another  problem  is  encountered  in  the  40  m  duct  height  case  when  a  large 
number  of  zeroes  are  found  in  the  lower  half  complex  qjj  plane.  These  zeroes 
appear  to  belong  to  the  reflection  coefficient  on  the  wrong  sheet  of  the  branch  cut 
and  are  not  waveguide  modes.  This  happens  because  the  search  region  has  been 
extended  below  the  real  q^  axis  to  avoid  the  singularity  in  SURF.  The  problem  with 
this  singularity  should  have  been  solved  within  SURF,  especially  because  it  occurs 
only  when  the  derivative  of  the  subroutine  output  variable  gamma  with  respect  to  qjj 
is  computed.  Since  this  derivative  is  not  needed  during  mode  search,  the  extension 
of  the  search  region  to  the  negative  q,,  plane  is  unnecessary.  A  simplified  routine, 
SURFO,  is  introduced  which  is  exactly  the  same  as  SURF  except  that  it  does  not 
evaluate  the  derivative  of  gamma.  By  using  this  subroutine  instead  of  SURF,  the 
search  path  in  the  revised  program  does  not  avoid  the  real  and  the  imaginary  axes. 

1.  FNDMOD 

The  search  region  is  limited  to  the  upper  half  qj,  plane.  All  the  modes 
found  are  ordered  according  to  their  range  attenuation  rates  before  those  numbered 
beyond  the  maximum  modes  allowed  are  abandoned. 

2.  FZEROX 

Since  the  curve  Im{F(q)}=0  enters  into  a  mesh  square  through  an  edge, 
the  values  of  Im{F(q)}  must  change  sign  over  the  end  points  of  either  one  or  all 
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three  other  edges.  When  there  is  only  one  other  edge  across  which  Im{F(q)}  changes 
sign  at  its  end  points,  it  is  the  edge  across  which  the  curve  Im{F(q)}=0  exits  the 
mesh  square.  Ambiguity  arises  when  all  edges  indicate  a  change  of  sign  at  their  end 
points.  When  this  occurs,  a  "right  turn  rule"  is  adopted  which  assumes  that  the  curve 
exits  the  edge  to  the  right  of  the  one  along  which  it  enters  the  mesh  square.  Such  a 
rule  avoids  the  retracing  of  the  search  path  when  the  mesh  square  is  revisited  as 
entering  this  same  mesh  square  from  the  left  side  of  an  edge  after  exiting  from  its 
right  side  requires  a  crossing  of  the  Im{F(q)}=0  curve,  which  is  prohibited  under  the 
simple  zero  assumption.  On  the  other  hand,  the  actual  curv'e  may  have  turned  left 
and  then  returns  to  this  mesh  square,  i.e., following  a  "left  turn  rule."  Under  such  a 
scenario,  this  wrong  choice  would  have  left  a  segment  of  the  curve  not  searched.  This 
difficulty  has  not  been  obsers'ed  during  testing.  In  fact  the  ambiguous  situation 
seldom  occurs.  Note  also  that,  as  remarked  above,  two  lines  of  Im{F(q)}=0  do  not 
cross  each  other  unless  a  higher  order  zero  is  present.  Hence  only  a  right  turn  rule 
or  a  left  turn  rule  for  the  curve  to  exit  the  mesh  is  allowed.  Exiting  the  opposite  edge 
demands  a  pair  of  crossing  Im{F(q)}=0  curves  within  the  mesh  square.  This  violates 
the  assumption  that  all  zeroes  are  simple.  Also  note  that,  the  possibility  of  vanishing 
Re{F(q)}  or  Im{F(q)}  values  at  the  comers  of  a  mesh  square  is  eliminated  through 
a  small  adjustment  in  FINDFX. 

3.  FINDFX 

Both  the  vertical  shift  away  from  the  real  q]]  axis  and  the  horizontal 
offset  away  from  the  imaginary  axis  are  unnecessary'  and  have  been  removed  from 
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this  routine.  Furthermore,  as  a  result  of  converting  to  the  complex  exponent 
representation,  the  sine  and  cosine  of  the  argument  of  the  modal  function  are 
examined  for  sign  changes  in  FZEROX.  This  is  implemented  in  FINDFX  by 
including  the  cosine  and  sine  values  of  the  argument  of  the  modal  function  in  the 
output  list.  To  avoid  the  indeterminate  case  when  either  the  real  or  the  imaginary 
part  of  the  modal  function  becomes  zero  at  any  comer  of  a  mesh  square,  the 
argument  for  computing  the  cosine  and  sine  values  is  increased  by  2“^^  w'hen  this 
occurs.  This  is  equivalent  to  a  consistent  small  distortion  of  the  particular  comer  of 
the  mesh  square.  This  will  not  cause  any  error  in  locating  the  zero  because  FINDFX 
still  returns  separately  the  unmodified  exponent  of  the  value  of  the  modal  function. 

4.  ROOTS 

Assuming  that  the  modal  function  is  analytic  within  the  mesh,  this 
subroutine  utilizes  the  values  of  the  modal  function  at  the  four  corners  of  the  mesh 
square  to  determine  the  Taylor  series  expansion  coefficients  of  the  modal  function 
to  the  third  order.  The  roots  of  this  cubic  polynomial  are  then  located  using  Cardan's 
solution  by  radicals.  If  the  higher  order  coefficients  fall  below  machine  resolution  for 
a  root  w'ithin  the  mesh  square,  these  coefficients  are  regarded  as  zero  and  the  order 
of  the  polynomial  is  reduced  and  can  be  solved  more  expediently.  If  the  function  is 
determined  to  be  constant  over  the  mesh  square,  the  center  of  the  square  is  taken 
as  the  root  location. 
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D.  EVALUATING/4,  AND 

As  discussed  in  the  Introduction,  the  /I,  and  B,  coefficients  can  be  evaluated 
either  from  the  top  level  down  or  from  the  lowest  level  up.  These  two  procedures  are 
simply  called  "integration  down"  and  "integration  up"  respectively  in  the  original 
documentation  [Ref.  4].  The  location  of  a  mode  has  been  called  an  eigenvalue.  That 
the  results  of  integration  down  and  integration  up  agree  is  a  manifestation  that  the 
eigenvalue  is  located  accurately. 

The  subroutine  ABCOEF  evaluates  the  coefficients  and  B,  for  each  mode. 
If  the  range  attenuation  rate  for  a  mode  is  greater  than  0.1  dB/km,  the  coefficients 
are  evaluated  from  the  lowest  layer  up.  Otherwise,  it  is  evaluated  tiom  the  top  layer 
down.  It  is  obvious  that  such  a  rule  must  be  implemented  because  the  results  of 
integration  up  and  integration  down  do  not  agree  for  many  modes.  Efforts  are  made 
to  determine  the  cause  of  this  discrepancy  and  to  devise  a  means  to  resolve  it. 

Investigation  reveals  that  inadequate  precision  in  the  location  of  the  modes  is 
one  source  of  the  problem.  Since  the  B,  coefficients  depend  on  the  coefficients 
while  the  A^  coefficients  are  obtained  directly,  only  the  A^  coefficients  need  to  be 
examined.  The  /4,  coefficients  of  the  six  modes  of  lowest  range  attenuation  rates  for 
all  21  profiles  except  the  one  without  evaporation  duct  are  computed  using 
eigenvalues  of  different  accuracy  controlled  by  the  first  order  Newton-Raphson 
iteration  method.  Table  1  shows  the  A^  coefficient  computed  with  the  new  program. 
They  are  arranged  from  the  top  layer  down.  In  the  i-th  layer,  the  A^  coefficient 
computed  by  integration  downward  depends  only  on  A^^j  in  the  layer  above  while 
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that  computed  by  integration  upward  depends  only  on  in  the  layer  below.  Hence 
in  each  layer,  the  coefficient  obtained  by  integration  downwaro  is  listed  above  that 
obtained  by  integration  upward.  There  are  five  sets  of  Aj  values  listed,  with  the 
magnitudes  given  in  powers  of  10  and  the  phase  given  as  a  multiple  of  r.  They  are 
obtained  from  eigenvalues  of  decreasing  accuracy,  the  one  used  to  compute  the  left 
most  column  being  the  most  accurate.  The  first  set  is  computed  using  an  eigenvalue 
having  a  relative  accuracy  of  2“'^®;  The  second  set  uses  an  eigenvalue  with  a  relative 
accuracy  of  2“^^;  The  relative  accuracy  of  the  eigenvalue  for  the  third  set  is  2“^^; 
For  the  fourth  set,  the  first  order  Newton-Raphson  iteration  of  the  mode  location  is 
set  at  an  absolute  accuracy  of  0.03  of  the  mesh  size,  sanv.  as  that  specified  in  the 
original  program;  The  eigenvalue  for  the  right  most  set  is  the  mode  location 
estimated  by  ROOTS  without  modification  by  the  Newton-Raphson  iteration.  It  is 
clear  that,  for  this  mode,  the  difference  between  these  two  methods  of  computing  the 
coefficients  becomes  negligible  as  the  accuracy  in  mode  location  increases.  For 
example,  in  the  8-th  layer,  the  magnitude  of  Aj  computed  by  integrating  downward 
changes  from  —1.9078  to  0.3482  to  0.3460  to  0.3459,  which  agrees  with  the  result 
computed  by  integrating  upward.  The  phase  follows  the  same  trend  to  an  agreement 
within  0.001  TT.  Table  2  shows  a  similar  set  of  output,  but  the  coefficients  fail  to  agree 
even  when  the  relative  accuracy  is  increased  to  Note  that  the  actual  difference 
in  both  the  real  part  and  the  imaginary  part  of  the  two  most  accurate  eigenvalues  is 
about  2““^*.  Double  precision  accuracy  appears  to  be  insufficient  for  the  coefficients 
computed  with  these  two  methods  to  agree  for  all  modes.  Some  interesting  features 
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can  be  observ'cd  in  both  tables,  which  are  present  in  all  120  sets  of  values  computed. 
When  disagreement  is  present  in  one  set  of  A-  coefficients  such  as  those  in  either 
Table  I  or  Table  2,  the  change  toward  smaller  differences  with  improving  eigenvalue 
accuracy  occurs  mainly  in  one  way  of  computation,  but  not  both.  For  example,  in 
Table  l,the  values  of  integration  downward  improve  with  better  eigenvalue  accuracy, 
while  those  computed  by  integrating  upward  change  little.  In  Table  2,  the  results  of 
integration  downward  are  the  ones  that  are  holding  steady  as  the  accuracy  in 
eigenvalue  improves.  Furthermore,  when  disagreement  occurs,  the  layer  in  which  the 
Aj  coefficient  has  the  smallest  magnitude,  i.e.,the  one  having  the  most  negative 
fxiwer  of  10,  divides  the  table  into  two  parts.  The  results  of  two  different  ways  of 
computation  agree  in  the  layers  above  this  one  if  they  disagree  in  those  below  it,  and 
vise  versa.  No  explanation  will  be  attempted.  Instead,  practical  rules  are  drawn  up 
to  take  advantage  of  these  facts.  In  Table  l,the  process  of  integration  upward  goes 
through  the  troublesome  10-th  layer  and  produces  results  which  agree  with  the  results 
of  downward  integration  before  the  downward  process  goes  through  the  10-th  layer. 
On  the  other  hand,  the  downward  integration  is  tripped  up  going  across  the  10-th 
layer  and  produces  results  which  fail  to  agree  with  the  results  from  upward 
integration.  It  is  clear  that  the  results  from  upward  integration  are  the  correct  ones. 
This  conclusion  is  further  supported  by  the  fact  that  improving  the  accuracy  of  the 
eigenvalue  does  not  change  significantly  the  results  of  upward  integration.  Similar 
argument  leads  to  the  conclusion  that  in  Table  2,  the  results  of  downward  integration 
are  the  correct  values. 
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It  can  be  concluded  from  the  above  observ'ations  that  one  of  the  methods  of 
computing  the  /I,  coefficients  converges  to  the  correct  value  much  faster  then  the 
other.  It  is  also  found  that  this  method  of  faster  convergence  is  always  able  to  arrive 
at  the  correct  values  for  /I,  for  all  the  cases  under  investigation. 

Table  3  lists  the  statistics  of  the  method  of  integration  which  yields  the  correct 
Ai  coefficients  for  each  of  the  120  modes  investigated.  The  differences  in  magnitudes 
and  phases  in  the  lowest  layer  and  in  the  layer  below  the  highest  are  also  listed. 
Since  for  most  of  the  cases  when  disagreement  in  A^  values  occurs,  the  correct 
integration  is  upward,  this  is  used  as  the  default.  To  decide  that  downward 
integration  should  be  utilized,  the  following  steps  are  taken:  The  first  A-  value  of 
downward  integration  is  computed  and  compared  to  the  value  from  upward 
integration.  If  the  magnitudes  in  dB  disagree  by  less  than  0.02  dB,  their  phases  will 
be  checked.  If  the  phases  differ  by  less  than  IO'^tt,  the  agreement  is  deemed 
acceptable  and  the  /I,  and  B,  coefficients  computed  from  the  low'est  layer  up  are 
used.  Otherwise,  the  coefficients  are  re-evaluated  again  from  the  highest  layer  down. 

Once  the  correct  method  of  evaluating  the  A-  and  coefficients  is  used,  the 
accuracy  of  the  mode  location  becomes  less  critical.  For  all  the  cases  investigated, 
the  Ai  coefficients  obtained  from  mode  locations  estimated  w'ith  or  without  the 
Newton-Raphson  first  order  iteration  differ  only  by  0.06  dB  in  magnitude  and 
O.OOlSrin  phase  at  most.  In  fact,  few  cases  show  differences  more  than  0.002  dB  and 
0.0001  X.  The  Newton-Raphson  iteration  is  not  needed.  Hence  the  subroutines 
NOMSHX,  FDFDTX  and  DXDETR  are  removed. 
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TABLE  3.  CONTINXED  2 
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TABLES.  CONTINUED  4 
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III.  CONCLUSION  AND  RECOMMENDATION 


A.  Performance 

This  revision  of  M-Layer  converts  the  extended  complex  number  representation 
of  an  exponentially  large  or  small  number  into  the  direct  representation  by  its 
complex  exponent.  The  accuracy  of  the  computation  has  been  improved  in  two  ways: 
First,  an  interpolation  algorithm  has  been  devised  when  severe  cancellation  of  the 
addends  is  detected.  Secondly,  accuracy  for  the  evaluation  of  the  Airy  function  has 
been  improved,  not  just  by  summing  the  Taylor  series  to  double  precision  resolution 
and  by  adopting  six-term  Gaussian  quadrature,  but  also  by  expanding  the  region 
within  which  the  more  expedient  Gaussian  quadrature  is  excluded  in  favor  of  the 
more  accurate  but  time-consuming  Taylor  series  summation.  The  improvement  in 
accuracy  is  most  easily  seen  from  Table  1. 

As  discussed  in  the  Introduction,  evaluating  the  A-  and  fi,  coefficients  either 
from  the  lowest  layer  up  (integration  up)  or  from  the  top  layer  down  (integration 
down)  must  result  in  the  same  values.  This  property  provides  a  consistency  check  for 
the  accuracy  of  the  computation.  For  the  six  modes  of  lowest  range  attenuation  rates 
of  the  20  profiles  of  different  duct  heights.  Table  1  lists  the  maximum  difference  for 
each  mode  which  shows  a  discrepancy  between  these  two  methods  of  evaluating  the 
/I,  coefficients.  For  each  profile,  the  maximum  value  in  magnitude  difference  in  dB 
among  all  the  layers  is  listed  if  it  is  greater  than  2.  If  the  phases  of  the  coefficients 


TABLE  1.  MAXIMUM  DIFFERENCE  IN  A,  COEFFICIENT  BETWEEN 
INTEGRATION  UP  AND  DOWN 


TABLE  1.  CONTINUED 


Duct 

height 

(m) 

Mode 

Difference  in  A,  coefficient 

Magnitude  difference  in 
(dB) 

Phase  difference  over 

0.1  T 

original 

revised 

original 

revised 

34 

2 

37.98 

Yes 

3 

715.7 

209.94 

Yes 

Yes 

36 

2 

112.74 

Yes 

3 

957.92 

231.68 

Yes 

Yes 

38 

2 

107.44 

52.26 

Yes 

Yes 

4 

1249 

255.8 

Yes 

Yes 

40 

3 

167 

112.72 

Yes 

Yes 

4 

823.56 

258.18 

Yes 

Yes 

Magnitude  difference  within  2dB  are  not  listed. 

deviate  more  than  O.lirin  any  layer,  that  particular  mode  is  also  singled  out.  The 
location  of  the  mode  of  the  revised  program  is  within  a  relative  accuracy  of  2“*^® 
achieved  through  first  order  Newton-Raphson  iteration.  Even  though  discrepancies 
still  exist  when  the  duct  is  28  meters  or  higher,  it  is  clear  that  the  revised  program 
computes  more  accurately  than  the  original  one. 

For  the  cases  where  the  two  methods  of  evaluating  the  and  B,  coefficients 
disagree,  it  has  been  observed  that  one  of  the  methods  always  leads  to  A-  values 
which  are  little  changed  when  the  accuracy  in  mode  location  is  varied,  while  the 
other  method  produces  /I,  values  which  shift  toward  the  results  of  the  other  method 
as  the  accuracy  of  mode  location  improves.  Based  on  this  observation,  a  consistency 


check  is  implemented  into  the  program  to  identify  the  method  which  converges 
better.  For  the  120  cases  investigated,  when  this  method  of  faster  convergence  is 
used,  the  coefficients  obtained  from  mode  locations  estimated  with  or  without  the 
Newton-Raphson  first  order  iteration  differ  only  by  0.06  dB  in  magnitude  and 
0.0013xin  phase  at  most.  In  fact,  few  cases  show  differences  more  than  0.002dB  and 
0.000 lx. This  allowed  the  Newton-Raphson  iteration  to  be  removed  in  this  revision. 

Table  2  compares  the  performance  between  the  original  and  the  revised 
programs.  The  time  spent  to  find  the  modes  has  been  reduced  by  an  average  of 
22.58%.  The  revised  program  can  always  produce  the  modes  found  by  the  original 
program.  Moreover,  the  mode  search  is  stable  for  the  new  program:  the  time  it 
requires  to  search  for  the  modes  is  about  the  same  for  similar  profiles.  The  sudden 
jumps  in  mode  search  time  for  the  24  m  and  the  40  m  cases,  which  indicate  troubles 
during  the  search,  no  longer  happen. 

With  the  proper  method  of  evaluating  the  and  B,  coefficients  determined  by 
the  consistency  check,  the  output  of  the  revised  program  differs  from  the  original 
program  in  some  cases.  The  most  serious  deviation  has  been  obser\'ed  for  the  38  m 
duct  height  case  as  shown  in  Tables  3  and  4.  For  example,  at  a  range  of  36.5  km  with 
the  transmitter  at  a  height  of  25  m  and  the  receiver  at  10  m,  the  coherent  path  loss 
is  175.93  dB  from  the  original  program,  and  is  167.90dB  from  the  revised  program. 
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TABLE  2.  OVERALLMODE  SEARCH  PERFORMANCE  COMPARISON 


DUCT 

HEIGHT 

(meters) 

ORIGINAL  PROGRAM 

REVISED  PROGRAM 

Time 

Improvement 

Time 

Modes 

Time 

Modes 

00 

0:00:37 

3 

0:00:35 

3 

5.40% 

02 

0:32:14 

9 

0:31:55 

9 

0.98% 

04 

1:14:12 

25 

1:05:04 

25 

12.31% 

06 

2:10:18 

53 

1:56:50 

53 

10.33% 

08 

0:35:58 

39 

0:29:25 

39 

18.21% 

10 

0:53:24 

59 

0:48:32 

61 

9.11% 

12 

1:09:40 

86 

1:01:44 

89 

11.39% 

14 

1:20:42 

94 

1:11:13 

97 

11.75% 

16 

1:54:35 

95 

1:18:07 

97 

31.82% 

18 

1:45:09 

100 

1:27:15 

104 

17.02% 

20 

1:46:19 

103 

1:34:20 

105 

11.27% 

on 

1:52:54 

105 

1:35:18 

lOo 

15.59% 

24 

3:42:59 

106 

1:46:47 

107 

52.11% 

26 

2:07:42 

106 

1:43:55 

108 

18.62% 

28 

2:00:05 

107 

1:44:59 

109 

12.57% 

30 

1:59:59 

107 

1:46:19 

108 

11.39% 

32 

1:55:29 

108 

1:42:58 

110 

10.84% 

34 

2:29:57 

109 

2:15:58 

111 

9.32% 

36 

2:31:40 

109 

2:17:20 

112 

9.45% 

38 

2:38:44 

110 

2:18:09 

111 

12.97% 

40 

5:41:17 

95 

2:39:39 

111 

53.22% 

Total 

40:23:54 

31:16:22 

22.58% 

TABLE  3.  ORIGINAL  PROGRAM  OUTPUT:  38  M  DUCT 


frequency  = 

9600. OOOC 

mhz 

range 

zt 

zr 

coherent 

incoherent 

coherent 

incoherent 

horizon 

(km) 

(m) 

(m) 

mode  sun 
(db) 

mode  sun 
(db) 

path  loss 
(db) 

path  loss 
(db) 

(km) 

27.3 

25.0 

4.0 

-15.30 

-15.62 

156.10 

156.43 

28.9 

27.3 

25.0 

6.0 

.62 

-2.35 

140.18 

143.16 

30.7 

27.3 

25.0 

8.0 

-1.11 

-4.21 

141.92 

145.01 

32.3 

27.3 

25.0 

10.0 

-27.26 

-12.66 

168.06 

153.46 

33.6 

36.5 

25.0 

4.0 

-16.94 

-16.62 

160.28 

159.96 

28.9 

36.5 

25.0 

6.0 

-.73 

-2.05 

144.07 

145.39 

30.7 

36.5 

25.0 

8.0 

-2.21 

-3.72 

145.55 

147.06 

32.3 

36.5 

25.0 

10.0 

-32.59 

-14.29 

175.93 

157.64 

33.6 

45.8 

25.0 

4.0 

-19.89 

-16.96 

165.20 

162.26 

28.9 

45.8 

25.0 

6.0 

-2.81 

-1.89 

148.11 

147.19 

30.7 

45.8 

25.0 

8.0 

-4.11 

-3.43 

149.41 

148.74 

32.3 

45.8 

25.0 

10.0 

-28.57 

-15.22 

173.88 

160.52 

33.6 

TABLE  4.  REVISED  PROGRAM  OUTPUT:  38  M  DUCT 


I  frequency  = 

9600. OOOC 

1  mhz 

range 

zt 

zr 

coherent 

incoherent 

coherent 

incoherent 

horizon 

(km) 

(m) 

(m) 

mode  sun 
(db) 

mode  sun 
(db) 

path  loss 
(db) 

path  loss 
(db) 

(km) 

27.3 

25.0 

4.0 

-14.38 

-15.66 

155.18 

156.47 

28.9 

27.3 

25.0 

6.0 

.42 

-2.37 

140.39 

143.18 

30.7 

27.3 

25.0 

8.0 

-1.52 

-4.21 

142.33 

145.02 

32.3 

27.3 

25.0 

10.0 

-21.20 

-12.51 

162.01 

153.31 

33.6 

36.5 

25.0 

4.0 

-17.32 

-16.60 

160.66 

159.94 

28.9 

36.5 

25.0 

-■.0 

-.48 

-2.08 

143.82 

145.42 

30.7 

36.5 

25.0 

8.0 

-1.62 

-3.73 

144.96 

147.07 

32.3 

36.5 

25.0 

10.0 

-24.56 

-14.04 

167.90 

157.38 

33.6 

45.8 

25.0 

4.0 

-20.26 

-16.93 

165.57 

162.23 

28.9 

45.8 

25.0 

6.0 

-3.14 

-1.93 

148.44 

147.23 

30.7 

45.8 

25.0 

8.0 

-4.62 

-3.46 

149.92 

148.76 

32.3 

45.8 

25.0 

10.0 

-25.40 

-14.90 

170.71 

160.21 

33.6 
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B.  Recommendation 


The  mode  search  protocol  of  this  program  needs  to  be  revised.  Since  the  search 
is  limited  by  the  maximum  range  attenuation  rate  accepted,  it  is  logical  to  begin  with 
locating  the  mode  of  the  lowest  or  the  highest  attenuation,  then  proceed  to  look  for 
the  next  one  in  the  order  of  increasing  or  decreasing  attenuation  rate.  Furthermore, 
under  the  assumption  of  analyticity  over  the  search  region,  there  should  be  only  one 
connected  "phase  line"  of  vanishing  real  part  of  the  modal  function  on  which  all  the 
modes  are  located.  The  partition  of  the  search  region  into  rectangles  as  has  been 
done  in  this  program  tends  to  cut  the  "phase  line"  into  segments  before  the  program 
starts  to  search  for  the  end  points  of  these  segments  and  then  follow  the  segments 
in  different  directions.  It  is  clear  that  a  better  way  is  to  search  for  one  end  of  the 
"phase  line"  along  a  line  of  a  constant  attenuation  rate  in  the  search  region,  either 
at  the  maximum  accepted  or  the  minimum  possible  attenuation,  then  follow  this 
"phase  line"  all  the  way  to  the  other  end.  This  technique  works  even  if  the  "phase 
line"  branches  off  into  several  directions  at  a  Stokes’  point. 
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APPENDIX  A:  SUBROUTINE  XCADD 


This  Appendix  lists  the  addition  subroutine  XCADD  which  returns  the  complex 
exponent  of  the  sum  when  the  complex  exponents  of  the  addends  are  given.  This  is 
a  complete  re-write  of  the  original  subroutine  of  the  same  name. 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 


subroutine  xcadd( zx, zlx, z2x) 
c 

c  Given  zlx  and  z2x,  this  subroutine  adds  the  two  complex  nunbers 
c  z1=exp(z1x)  and  z2=exp(z2x)  for  z=exp(zx)  and  returns  zx. 

c 

c  inputs... 

c  z1x=complex  exponent  of  the  complex  nur;.^er  z1 

c  z2x=coniplex  exponent  of  the  complex  nunber  z2 

c 

c  outputs... 

c  zx=complex  exponent  of  the  complex  nunber  z 

c 

c  subroutines  called... 

c 

c********** 

implicit  real*8  (a-h,o-z) 

comp I  ex* 16  zx, zlx, z2x, z tlx, zt2x, clo9zh,dsum,c zero, c err x, cone, chpi 
parameter (pi =3. 141592653589793238462643d0, twopi =2.d0*pi , 

+  hpi =0.5d0*pi , zero=0.d0,c16=1 .d0/6.d0, 

♦  bit14=1 .d0/16384.d0,bit24=bit14/1024.d0,ctol=bit14, 

+  dpi =2259. d0/4294967296.d0/4294 967296. d0,hdpi=dpi/2.d0, 

+  e2m54=-3.742994775023704819d1,e2p27=-0.5d0*e2m54, 

+  Chpi=(0.d0,1.57079632679489661923132d0),cone=(1.d0,0.d0), 

♦  czero=(0.d0,0.d0),cerrx=(-3.742994775023704819d1,0.d0)) 

c  cerrx=e2tn54=-54*log(2)=exponent  below  machine  accuracy 

dimension  ztmp(2) ,stmp(2) 
equivalence  (ztmp,clogzh),(stmp,dsum) 

^***** 

c  Replace  the  input  variables  with  a  local  variable  so  that 
c  equations  in  the  form  of  y=x+y  will  not  lead  to  confusion, 
c 

zt1x=z1x 

zt2x=z2x 

c 

clogzh=0.5d0*(zt1x-zt2x) 
dxh=ztmp( 1 ) 

if(dxh  .It.  zero)  then 
zx=zt2x 
dxh=-dxh 
else 

zx=zt1x 
end  i  f 
c******** 

c  machine  accuracy  =  2**(-53) 

c  2**(27)=e**e2p27 

c 

if  (dxh  .ge.  e2p27)  then 
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A8 

49 

50 

51 

52 

53 

54 

55 

56  c 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67  c 

68  c 

69  c 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79  c 

80 


return 

else 

zx=0.5d0*(zt1x+zt2x) 
dsijii=cdexp(  cl  ogzh ) 
dSLin=1  .dO/dsun+dsum 
if  (cdabsCdsLin)  .gt.  ctol)  then 
zx=cdlog(dsutn)+zx 
else 

Cancellation  is  serious.  Imtclogzh]  is  close  to  pi/2  or  -pi/2. 
yi=dnint(ztmp(2)/twopi )*2.d0 
ztmp(2)=ztmp(2)-pi*yi 
dyi =dpi*yi 

if  (ztmp(2)  .It.  zero)  then 
clogzh=-clogzh 
dyi=-dyi 
end  i  f 

ztmp<2)=(ztmp(2)-hpi )-hdpi'dyi 
dsui>=2.d0*c  logzh*(cone*c16*clogzh*clogzh) 
if  (dsum  .eq.  czero)  then 

Note  that  a  complete  cancellation  of  two  nonzero  numbers  of 
order  one  is  considered  to  be  as  accurate  as  what  is  allowed 
by  the  machine  and  the  algorithm. 
zx=cerrx»chpi+zx 
else 

dSLin^cd  log  (dsum) 

if  (stiipd)  .It.  e2m54)  stmpO )=e2m54 
zx=dsum+chpi+zx 
end  i  f 
end  i  f 
return 
end  i  f 

end 
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APPENDIX  B:  SUBROUTINE  FZEROX 


This  Appendix  includes  the  listing  of  the  subroutine  FZEROX  which  identifies 
the  meshes  which  may  contain  modes  within  a  contour  rectangle.  The  Shellman- 
Morffit  mode  locating  algorithm  has  been  completely  replaced. 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

n 

12 

13 

u 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

1 

2 

3 

4 

5 

6 
7 

32 

33 

34 

35 

36 

37 

38 


subroutine  f zerox( t lef t , tright, tbot , ttop, tmshO, zeros, ni , nf ) 

c  fzerox  is  a  routine  for  finding  the  zeroes  of  a  complex  function,  f, 
c  which  lie  within  ^  specified  rectangular  region  of  the 

c  complex  q11  plane,  assuming  that  the  function  has  onl/ 

c  simple  zeroes  over  this  rectangle, 
c 

c  parameters 
c  tleft  - 

c  tright- 

c  tbot  - 

c 

c  ttop  - 

c  tmesh  - 

c 
c 
c 

c  zeros  - 

c 

c  nf-ni  - 

c 

c  subroutines  calledd-- 
c  findfx 

c  roots 

c  nomshx 

implicit  double  precision  (a-h,o-z) 

co(iplex»16  f10,f01,f11,fxnew,fxold,fxC0,fx10,fx01,fx11, 

♦  czero, one,c i , sol , zeros 

parameter(czero=(0.d0,0.d0) , one=( 1 .dO, O.dO) ,c i =(0 .dO, 1 .dO) ) 
Sinclude:  'mlaparm. inc ' 

'***  Begin  listing  of:  mlaparm. inc 
c 

c  include  file  to  define  the 

c  maximur  M  of  layers  (mxlayr) 

c  maximum  #  of  modes  (mxmode) 

c 

parameter  (mxlayr=35  ) 
parameter  (mxmode=127) 

'***  End  listing  of:  mlaparm. inc 

dimens i on  kedgel < 100) , kedge2( 100) , kedge3( 100) , kedge4( IOC) , 
c  ♦  loc12r(inxmode) ,  loc12i  (mxmode) ,  loc?3r(mxniode) ,  loc23  i  (raxmode) , 
c  ♦  loc34r(mxmode) , loc34i (mxmode) , loc41 r(mxmode) , loc41 i (mxmode) , 

♦  sol (3), theta(2) , zeros(2*mxmode*1 ) 
c 

c 

common  /tmccom/tmesh 


specifying  the  search  rectangle: 
value  of  the  real  part  of  q11  at  the  left  edge, 
value  of  the  real  part  of  q11  at  the  right  edge, 
value  of  the  imaginary  part  of  q11  at  the  bottom  edge, 
(this  is  set  to  0. ) 

value  of  the  imaginary  part  of  q11  at  the  top  edge, 
set  equal  to  about  half  the  average  spacing  between 
zeroes  within  the  rectangle.  A  smaller  value  may  be  used 
as  a  safety  measure,  but  tco  small  a  value  will  result 
in  excessively  long  run  time. 

output  list  ot  (complex)  values  of  q1 1  at  which 

zeroes  are  found. 

the  ncmber  of  zeroes  found 
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39  c***** 

40  c  maxnsq  -  maximun  number  of  mesn  squares  allowed  on  any  one 

41  c  phase  line 

42  c  maxnt  -  maximun  number  of  times  fzerox  will  reduce  tmesh 

43  c 

44  niaxnsq=3*max0(  int( (t top- tbot)/tmsh0) ,  int( ( tri ght-t  left  )/tmsh0) ) 

45  maxnt=2 

46  0***** 

47  tmesh  -  tmshO 

48  ntime  =  0 

49  go  to  7 

50  c 

51  5  tmesh=tmesh/2.0d0 

52  ntime  =  ntime+1 

53  ifCntime  .gt.  maxnt)  go  to  97 

54  c 

55  7  continue 

56 

57  c***** 

58  c  calculate  coordinates  of  rectangle  edges  in  tmesh  units 

59  c 

60  jit  =  idnintftlef t/tmesh-0.5d0) 

61  jrt  =  idni nt ( t r i gh t/tmesh*0.5d0) 

62  jtop  =  idn  int(ttop/ tmesh*  1 .5dC:) 

63  jbot  =  0 

64  c 

65  c  initialize  parameters  for  starting  search  at  upper  left 

66  c  corner  of  search  rectangle 

67  c 

68  ki  =  jtop 

69  k  r  =  j  1 1 

70  kedge  =  1 

71  call  f indf x( kr , k i , f xnew, xnew, ynew) 

72  nre1=0 

73  nre2=0 

74  nre3=0 

75  nre4=0 

76  knot12=0 

77  knot23=0 

78  knot34=0 

79  knot41=0 

SO  fir*ni 

81  ni1=ni*1 

82  go  to  15 

83  c . 

84  10  continue 

85  if(nrzl  .It.  2)  go  to  15 
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86 

c 

write(16,2000)  nrzl 

87 

go  to  5 

88 

15 

nrzl=0 

89 

nrsqu  =  0 

90 

20 

fxold=fxneH 

91 

xold=xnew 

92 

yold=yneu 

93 

go  to  (21 ,26,31 ,36) ,kedge 

94 

95 

c 

search  along  left  edge  of  rectangle  for  changes  tn  the 

96 

c 

sign  of  imag(f) 

97 

c 

98 

21 

continue 

99 

if (ki .eq. jbot)  then 

100 

kedge=2 

101 

go  to  26 

102 

end  i  f 

103 

ki  =  ki  - 1 

104 

cal  1  f indfx(kr,ki ,fxnew,xnew,yneH) 

105 

if  (yold*ynew  .gt.  O.dO)  go  to  20 

106 

i f(nre1 .eq.O)  go  to  23 

107 

c 

108 

c 

check  if  crossing  point  has  been  previously  found 

109 

c 

110 

do  22  k=1 ,nre1 

111 

if(ki .eq.kedgel(k))  go  to  20 

112 

22 

continue 

113 

c 

114 

c 

follow  phase  line  through  rectangular  region 

115 

c 

116 

23 

fxOI =fxold 

117 

f xOI r=xold 

118 

f xOI i =yold 

119 

f x00=f xnew 

120 

fx00r=xnew 

121 

fx00i=ynew 

122 

li  =  ki 

123 

Ir  =  jit 

124 

go  to  43 

125 

Qliimitit 

126 

C 

search  along  bottom  edge  of  rectangle  for  changes  in  the 

127 

r 

sign  of  imag(f> 

128 

C 

129 

26 

continue 

130 

if(kr.eq. Jrt)  then 

131 

kedge=3 

132 

go  to  31 
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133  end  if 

134  kr  =  kr+1 

135  call  f indfx(kr, ki , fxnew, xneu, ynew) 

136  if  (yold*yneu  .gt.  O.dO)  go  to  20 

137  i f (nre2.eq.0)  go  to  28 

138  c 

139  c  check  if  crossing  point  has  been  previously  found 

140  c 

141  do  27  k=1,nre2 

142  if(kr.eq.kedge2(k))  go  to  20 

143  27  continue 

144  c 

145  c  follow  phase  line  through  rectangular  region 

146  c 

147  28  fx00=fxold 

148  fx00r=xold 

149  fx00i=yold 

150  fx10=fxneu 

151  fx10r=xnew 

152  fx10i=ynew 

153  li  =  jbot 

154  Ir  =  kr-1 

155  go  to  48 

156  c***»» 

157  c  search  along  right  edge  of  rectangle  for  sign  changes  in  imag(f). 

158  c 

159  31  continue 

160  i f (ki .eq. j top)  then 

161  kedge=4 

162  go  to  36 

163  erxd  if 

164  ki  =  ki+1 

165  call  f indf x(kr , ki , f xnew, xnew.ynew) 

166  if  (yold*ynew  .gt.  O.dO)  go  to  20 

167  if(nre3.eq.O)  go  to  33 

168  c 

169  c  check  if  crossing  point  has  been  previously  found 

170  c 

171  do  32  k=1,nre3 

172  i f (ki .eq. kedge3{k ) )  go  to  20 

173  32  continue 

174  c 

175  c  follow  phase  line  through  rectangular  region 

176  c 

177  33  fx10=fxold 

178  fx10r=xold 

179  fxIOiryold 
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180  fx11=fxneH 

181  fx11r=xneK 

182  fx11i=ynew 

183  li  =  ki-1 

184  Ir  =  jrt-1 

185  go  to  53 

186 

187  c  search  along  top  edge  of  rectangle  for  sign  changes  in  imag(f). 

188  c 

189  36  continue 

190  ifCkr.eq. j It)  go  to  80 

191  kr  =  kr-1 

192  call  f indfxCkr, ki , fxnew,xnew,ynew) 

193  if  (yold*ynew  .gt.  O.dO)  go  to  20 

194  if(nre4.eq.0)  go  to  38 

195  c 

196  c  check  if  crossing  point  has  been  previously  found 

197  c 

198  do  37  k=1 ,nre4 

199  i f (kr .eq.kedge4( k ) )  go  to  20 

200  37  continue 

201  c 

202  c  follow  phase  line  through  rectangular  region 

203  c 

204  38  fxll^fxold 

205  fx11r=xold 

206  fx11i=yold 

207  fx01=fxnew 

208  fx01r=xnew 

209  fx01i=ynew 

210  li  =  jtop-1 

211  Ir  =  kr 

212  go  to  58 

213  c****» 

214  c  enter  mesh  square  from  left  side  or  exit  rectangle  at  right  edge. 

215 

216  41  lr=lr+1 

217  if  (Ir  .le.  jrt-1)  go  to  42 

218  nre3=nre3+1 

219  kedge3(nre3)=l i+1 

220  go  to  10 

221  42  fx01=fx11 

222  fx01r=fx11r 

223  fx01i=fx11i 

224  fx00=fx10 

225  fxOOrxfxIOr 

226  fx00i=fx10i 


46 


227 

43 

continue 

228 

cal  1  f indfxC lr+1 , 1 i  +  1 , fxl 1 , fxl 1 r, fxl  1 1 ) 

229 

call  f indfxC lr+1,li,fx1Q,fx10r,fx10i) 

230 

c** 

***** 

231 

c 

Determine  the  edge  of  exit  of  im(f)=0  from  current  mesh. 

232 

edgei t=fx01 i*fx1 1 i 

233 

edgeib=fx00i*fx10i 

234 

if  (edgeib  .gt.  O.dO)  then 

235 

c 

lm(f)=0  goes  through  the  01  to  10  line. 

236 

if  (edgei t  .gt.  O.dO)  then 

237 

c 

lm(f)=0  goes  through  the  10  to  11  edge  (edge  1). 

238 

lout=1 

239 

else 

240 

c 

lm(f)=0  goes  through  the  01  to  11  edge  (edge  2) 

241 

lout=2 

242 

end  i  f 

243 

else 

244 

c 

lm(f)=:0  goes  through  the  00  to  10  edge  (edge  4) 

245 

lout=4 

246 

if  (edgeit  .It.  O.dO)  then 

247 

c 

lm(f)=0  also  runs  through  01  to  11  and  10  to  11  edges. 

248 

c 

Store  crossing  location  and  in/out  information. 

249 

knot34=knot34+1 

250 

c 

loc34r(knot34)=lr 

251 

c 

loc34i (knot34)=l i 

252 

end  if 

253 

end  i  f 

254 

^*****«*  1 

255 

go  to  60  1 

256 

c*****  1 

257 

c 

enter  mesh  square  from  bottom  side  or  exit  rectangle  at  top  edge. 

258 

46 

li=li+1 

259 

if  (li  .le.  jtop-1)  go  to  47 

260 

nre4=nre4*1 

261 

kedge4(nre4)= 1 r 

2'  . 

go  to  10 

263 

47 

fx00=fx01 

264 

fx00r=fx01r 

265 

fx00i=fx01 i 

266 

fx10=fx11 

267 

fx10r=fx11r 

268 

fx10i=fx11i 

269 

48 

continue 

270 

cal  1  f indfx( Ir,li  +  1,fx01,fx01r,fx01i) 

271 

cal  1  f indfx( lr+1 , 1 i  +  1 , f x1 1 , fxl 1r, fxl 1 i ) 

272 

^*******  1 

273 

c 

Determine  the  edge  of  exit  of  im(f)=0  from  current  mesh. 

47 

274 

edgei I  =  fx00i*fx01  i 

275 

edgei r=fx10i*fx1 1  i 

276 

if  (edgei r  .gt.  O.dO)  then 

277 

c 

lm(f)=0  goes  through  the  00  to  1'  iine. 

278 

if  (edgei 1  .gt.  O.dO)  then 

279 

c 

lm(f)=0  goes  through  the  01  to  11  edge  (edge  2) 

280 

lout=2 

281 

else 

282 

c 

lm(f)=0  goes  through  the  00  to  01  edge  (edge  3). 

283 

lout =3 

284 

end  i  f 

285 

else 

286 

c 

lm(f)=0  goes  through  the  10  to  11  edge  (edge  1) 

287 

lout=1 

288 

if  (edgei 1  .It.  O.dO)  then 

289 

c 

lm(f)=0  also  runs  through  00  to  01  ana  01  to  11  edges. 

290 

c 

Store  crossing  location  and  in/out  information. 

291 

ltnot41=knot41  +  l 

292 

c 

loc41 r(knot41 )= 1 r 

293 

c 

loc41 i (knot41  )  =  l i 

294 

end  i  f 

295 

end  i  f 

296 

c** 

***** 

297 

go  to  60 

298 

c** 

*** 

299 

c 

enter  mesh  square  from  right  side  or  exit  rectangle  at  left  edge 

300 

301 

51 

lr=lr-1 

302 

if  ( 1 r  ,ge.  jit)  go  to  52 

303 

nrel =nre1+1 

304 

kedgel (nrel )=l i 

305 

go  to  10 

306 

52 

fx1'i  =  fx01 

307 

fx11r=fx01r 

308 

fx11i=fx01i 

309 

fx10=fx00 

310 

fx10r=:fx00r 

311 

fx10i=fx00i 

312 

53 

continue 

313 

call  findfx(lr,li+1,fx01,fx01r,fx01i) 

314 

call  findfx(lr,li, f xOO, f xOOr, f xOOi ) 

315 

^*«****« 

316 

c 

Determine  the  edge  of  exit  of  im(f)=0  from  current  mesh. 

317 

edgeit=f xOI i*fx1 1 i 

318 

edgei b=fx00i*fx10i 

319 

if  (edgeit  .gt.  O.dO)  then 

320 

c 

lm(f)=0  goes  through  the  01  to  10  line. 
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321 

322 

323 

324 

325 

326 

327 

328 

329 

330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 

361 

362 

363 

364 

365 

366 

367 


if  (edgeib  .gt.  O.dO)  then 

c  lm(f)=0  goes  through  the  00  to  01  edge  (edge  3). 

lout=3 
else 

c  lm(f)=0  goes  through  the  00  to  10  edge  (edge  4) 

lout=:4 
end  i  f 
else 

c  lm(f)=0  goes  through  the  01  to  11  edge  (edge  2) 

lout=2 

if  (edgeib  .It.  O.dO)  then 

c  lm(f)=0  also  runs  through  00  to  10  and  00  to  01  edges, 

c  Store  crossing  location  and  in/out  information. 

knot12=knot12+1 
c  loc12r(knot12)=lr 

c  loc12i (knot  12)  =  I i 

end  i  f 
end  i  f 
^****«** 

go  to  60 

c  enter  mesh  square  from  top  side  or  exit  rectangle  at  bottom  edge. 

56  I  i  =  I  i • 1 

i '  ( I i  -ge.  jbot )  go  to  57 
nre2=nre2+ 1 
kedge2(nre2)= I r*1 
go  to  10 

57  fx01=fx00 
fx01r=fx00r 
fxOI i=fx00i 
fx11=fx10 
fx11r=fx10r 
fxll i=fx10i 

58  continue 

call  findfx(lr,li,fx00,fx00r, f xOQ  i ) 
call  f indfx( I r*1 , I i , fxlO, fxIOr, fxIOi ) 

^******* 

c  Determine  the  edge  of  exit  of  im(f)=0  from  current  mesh, 
edgei I =f xOOi *f xOI i 
edgei r=fx10i*fx1 1  i 
if  (edgeil  .gt.  O.dO)  then 
c  lm(f)=0  goes  through  the  00  to  11  line, 

if  (edgeir  .gt.  O.dO)  then 

c  lm(f)=0  goes  through  the  00  to  10  edge  (edge  4) 

lout=4 
else 

c  lm(f)=0  goes  through  the  10  to  11  edge  (edge  1). 
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lout=1 


368 

369  end  i f 

370  else 

371  c  lm(f)=0  goes  through  the  00  to  01  edge  (edge  3) 

372  lout=3 

373  if  (edgeir  .It.  O.dO)  then 

374  c  lm(f)=0  also  runs  through  00  to  10  and  10  to  11  edges. 

375  c  Store  crossing  location  and  in/out  information. 

376  knot23=knot23+1 

377  c  loc23r(knot23)=lr 

378  c  loc23i (knot23)=l i 

379  end  i f 

380  end  i f 

381  c 

382  c******* 

383  60  continue 

384  nrsqu=nrsqu+1 

385  ifCnrsqu  .gt.  maxnsq)  go  to  95 

386  c****** 

387  c  Test  for  there  being  at  least  one  re(f)=:0  line  entering  and 

388  c  leaving  the  mesh  square. 

389  c 

390  if  ((fx00r*fx10r  .gt.  O.dO)  .and.  (fx01r*fx11r  .gt.  O.dO) 

391  ♦  .and.  (fx00r*fx01r  .gt.  O.dO))  go  to  (41,46,51,56)  lout 

392  c 

393  c  Coniputate  the  values  of  the  modal  function  at  the  corners  of  a 

394  c  a  mesh  square  to  determine  its  Taylor  series  to  the  3rd  order 

395  c  for  estimating  its  root  locations. 

396  c 

397  c  f00=one 

398  f 10=cdexp( f xIO-f xOO) -one 

399  f 01=cdexp(fx01 - f x00)-one 

400  f 1 1=cdexp(fx1 1  - f xOO) -one 

401  c 

A02  ^♦♦♦*************************************<Hk»****'**'*H-'*>HHHk1k**»r********** 

403  c  write  (16,3001)  ni , nf , I r , I i , knot 12,knot23, knot34 , knot41 

404  c  3001  format(/'  ni,  nf,  Ir,  li  and  knot12,  23,  34  and  43  before  ROOTS 

405  c  ♦  2i6,2x,2i6,2x,4i6) 

406  c 

407  c***********  estimate  locations  of  zeroes  by  radicals  **•*»**»***»***** 

408  c 

409  call  roots(f 10, fOI , f 1 1 ,sol ,nrsol ) 

410  c 

411  do  63  n=1 ,nrsol 

412  ureal  =  dreal (sol (n) ) 

413  uimag  =  dimag(sol (n) ) 

414  if  (ureal  .It.  O.dO  .or.  ureal  .gt.  I.OdO)  gc  to  63 
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415  if  (utmag  .It.  O.dO  .or.  uimag  .gt.  I.OdO)  go  to  63 

416  62  theta(1)=( Ir+ureal )*tmesh 

417  theta(2)=( I i+uimag)*tmesh 

418  nf  =  nf+1 

419  zerosCnf )=dcmplx(theta(1 ), theta(2)) 

420  nrzl=nrzl+1 

421  63  continue 

422  c****** ****************************************************** ******* 

423  c  unite  (16,3002)  ni,nf,nrsol 

424  c  3002  formate/'  out  of  ROOTS  at  63,  ni,  nf  and  #  of  roots  ',3i4) 

425  c***************************************«*«****************4******** 

426  c  continue  following  the  phase  line 

427  go  to  (41,46,51,56)  lout 

428  c****** 

429  cc 

430  80  continue 

431  c 

432  return 

433 

434  95  continue 

435  write(16,950C) 

436  write(16,4001)lr,li,ni,nf, tmesh 

437  write(*  ,9500) 

438  4001  formatCgo  to  5  from  95  at  Ir,  li  =  ’ ,  i6, ' , ' ,  i6, '  ni,  nf  =‘,i6, 

439  ♦'('.'6,',  mesh  size  =',d14.6) 

440  go  to  5 

441  c***** 

442  97  continue 

443  uri te(16,9700) 

444  ur i te( 16,4002) lr,li,ni,nf, tmesh 

445  write(*  ,9700) 

446  4002  format('go  to  5  from  97  at  Ir,  li  = ' , i6, ' , ' , i6, '  ni,  nf  =',i6, 

447  ■*■',','6,',  mesh  size  =' ,dl4.6,/ '  zeroes  found  are  kept.') 

448  c  nf=ni 

449  c 

450  return 

451  c 

452  c****  format  statements 

453  9500  format(/5x, ' too  many  squares  on  same  phase  line  --  ', 

454  $  'reduce  tmesh  and  start  over') 

455  9700  format(/5x, ' tmesh  has  been  reduced  but  problems  remain  in', 

456  $  '  executing  fzerox') 

457  c 

458  end 
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APPENDIX  C:  SUBROUTINE  ROOTS 


This  Appendix  contains  the  listing  of  the  subroutine  ROOTS.  This  subroutine 
replaces  the  portion  of  the  subroutine  FZEROX  where  the  coefficients  of  a  quadratic 
equation  are  determined,  and  the  subroutine  QUAD  for  locating  the  zeroes  of  a 
quadratic  polynomial.  In  the  revised  subroutine  FZEROX,  the  roots  of  a  cubic 
polynomial  has  to  be  found.  This  subroutine  determines  these  zeroes  by  radicals. 
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1 

2 

3 

A 

5 

6 

7 

8 
9 

10 

11 

12 

13 

U 

15 

16 

17 

18 

19 

20 
21 
22 
23 
2A 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 


subroutine  roots  ( f 1 , f2, f3,sol ,nrsol ) 

q*************************************************4***«**««-*«********** 

c  This  subroutine  finds  the  roots  of  a  third  order  polynomial  by 
c  radicals  when  the  values  of  this  polynomial  at  z=0,  z=1,  z=i  and 
c  z=1+i  are  given  as  f0=1,  fl+fO,  f2+f0  and  f3»f0  respectively, 
c  Note  that  this  algorithm  takes  cubic  roots  of  two  complex  numbers 
c  (hence  the  name  'solution  by  radicals')  and  use  their  linear 
c  combinations  as  the  roots  of  a  third  order  polynomial. 

Q  *******  Ik  ***************************************************■*********<*  * 

implicit  real*8  (a-h,  o-z) 

complex* 16  f 1 , f2, f3, zero,one,ci ,ep14,em14,ep23,em23, 

+  fa, fb, fc,fd, fal , fa2, fa3, fa1s,p,q,del t , z, zm, u, v, sol 

parameter  (xbi t52=52.d0*0.69314718055994531d0, thrd=1 .d0/3.d0, 

+  bit50=1.d0/33554432.d0/33554432.d0,bit51=bit50/2.d0, 

+  bit52=bit51/2.d0,tol=0.001d0, 

♦  zero=C0.d0,0.d0),one=(1 .d0,0.d0),ci=(0.d0, 1 .dO), 

♦  ep14=(0.5d0,0.5d0),em14=(0.5d0,-0.5d0), 

+  ep23= ( -  0 . 5d0 , 0 . 86602540378443864675d0 ) , 

+  em23= ( -  0 . 5d0 , -  0 . 866025403784438646  75d0 ) ) 

dimension  sol(*) 
fa=one 

fb=(f2-ci*f 1+em14*f3) 
fc  =  ( (ep14»one)*f 1  - (em14+one)*f2+ci*f3) 
fd=(em14*(f2-f1)-ep14*f3) 
if  (cdabs(fb)  .le.  bitSO)  fb=zero 

if  (cdabs(fc)  .le.  bit51)  fc=zero 

if  (cdabs(fd)  .le.  bit52)  fd=zero 

if  (fd  .ne.  zero)  then 
f a1 =( - thrd)*f c/fd 
f a2=f b/fd 
f a3=f a/fd 
fa1s=fa1*fa1 
p=thrd*f a2- f a1 s 

q=0.5d0*(fa3+fa1*fa2)-fa1*fa1s 
if  (p  .eq.  zero)  then 
if  (q.  eq.  zero)  then 
nrsol=l 
sol ( 1 )=f a1 
return 
else 

nrsol=3 

u=<{-2.d0)*q)**thrd 
sol ( 1 )=u+fa1 
sol(2)=ep23*u+fa1 
sol(3)=em23*u*fa1 
return 
end  i  f 
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48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71  100 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 


else 

if  (q.  eq.  zero)  then 
nrsot=3 
sold  )=fa1 
u=cdsqrt((-3.d0)*p) 
sol(2)=fa1+u 
sol<3)=fa1 -u 
return 
else 
v=p/q 
z=p*v*v 
absz=cdabs(z) 
if  (absz  .It.  tol)  then 
zm=-z 

fn=dint(1 .d0-xbit52/dlog(absz)) 

lastn^idint ( fn)- 1 

dnn=fn-0.5d0 

dnd=fn'»1  .OdO 

del t=one 

do  100  nt=1,lastn 
dnn=dnn- 1 .dO 
dnd=dnd- 1 .dO 

del  t  =  (dnn/dnd)*del  t*ziiHone 
continue 

del t=(0.5d0*del t/q)**thrd 
u=p*delt 
v=-1 .dO/delt 
else 

de 1 1 =cdsqr  t ( one+  z ) -  one 
u=(q*del t )**thrd 
v=-p/u 
end  i  f 
nrsol=3 

sol ( 1 )=u+v*f a1 
sol (2)=ep23*u*em23*v+fa' 
sol (3)=em23*u*ep23*v*f al 
return 
end  i  f 
end  i  f 

else  if  (fc  .ne.  zero)  then 
if  (fb  .eq.  zero)  then 
if  (fa  .eq.  zero)  then 
nrsol=1 
sol ( 1 )=zero 
return 
else 

nrsol=2 
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95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115  200 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 


z=cdsqrt( - fa/f c) 
soUl  )  =  z 
soU2)  =  -z 
return 
end  i  f 
else 

fa1=0.5d0»fb/fc 

fa2=fa/fc 

z=fa2/fa1/fa1 

absz=cdabs(z) 

if  (absz  .It.  tol)  then 

fn=dint(1 .d0-xbit52/dlog(absz)) 

lastn=idint(fn)- 1 

dnn=fn-0.5d0 

dnd=fn+1 .OdO 

del t=one 

do  200  nt=1 , lastn 
dnn=dnn- 1 .dO 
dnd=dnd- 1 .dO 

del  t  =  (dnn/dnd)’del  t*2'*one 
continue 

delt=-0.5d0*delt/fal 

nrsol=2 

sol<1)=fa2*delt 
sol(2)=1 .dO/delt 
return 
else 

del t=cdsqrt(one- z) 
nrsol=2 

sol(1)=-fa1*(one-delt) 
sol(2)=-fal*<one*delt) 
return 
end  i  f 
end  i  f 

else  if  (fb  .ne.  zero)  then 
nrsol=1 
sol ( 1 )=- f 3/f b 
return 
else 


nrsol=1 
sol ( 1 )=ep14 
return 


APPENDIX  D:  SUBROUTINE  ABCOEF 


This  Appendix  contains  the  listing  of  the  subroutine  ABCOEF.  The  consistency 
self-checking  procedure  has  been  implemented  to  determine  the  correct  method  to 
evaluate  the  and  B,  coefficients. 
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subroutine  abcoef (zero.m) 


1 

2  ^***** 

3  c  For  each  mode  m,  this  suboutine  calculates  A-B  coefficients  in 

4  c  all  layers  for  combining  two  linearly  independent  solutions  of 

5  c  Stokes'  equation  to  form  the  height  gain  function: 

6  c 

7  c  height  gain=expCbcoefx( t,m))*<k1*exp(acoefx( l,m))»k2) 

8  c 

9  c  where  k1  and  k2  are  two  independent  solutions  to  Stokes' 

10  c  equation.  In  the  top  layer  (i.e.  nzlayr)  the  height  gain  is: 

11  c 

12  c  height  gain=exp{bcoef x( I ,m))*h2 

13 

14  c  where  h2  is  a  solution  to  the  Stokes'  equation  associated 

15  c  with  outgoing  energy  flow.  Here  k1  and  k2  are  proportional 

16  c  to  the  k1  and  k2  used  by  Marcus  and  the  h2  is  proportional 

17  c  to  a  modified  Hankie  function  of  order  1/3. 

18 

19  c  inputs... 

20  c  zero-an  eigenvalue  in  q11  space 

21 

22  c  outputs. . . 

23  c  acoefx-two  dimensional  array  of  complex  exponents 

24  c  coefficients  used  to  combine  two  linearly 

25  c  independent  solutions  of  stokes'  equation 

26  c  bcoefx-two  dimensional  array  of  complex  exponents 

27  c  coefficients  used  for  normalizing  the  height  gains 

28 

29  c  note:  acoefx  and  bcoefx  are  passed  by  the 

30  c  conmon  block  /pap2/ 

31 

32  c  subroutines  called... 

33  c  xcdai 

34  c  xcadd 

35 

36  c  conrran  block  areas... 

37  c  comi 

38  c  co(ii2 

39  c  papi 

40  c  pap2 

41  c***** 

42 

43  implicit  real*8(a-h,o-z) 

44  complex*16  acocfx, bcoefx, cqi J ,h2xq1 ,dh2xq1 ,h2xq2,dh2xq2, k1 xql , 

45  *  dklxql ,k1xq2,dk1xq2, k2xq1 ,dk2xq1 , k2xq2,dk2xq2,  h2dk  1 X, 

46  t  C#i2k1x,h2dk2x,dh2k2x,n(jnax,dcnax,numbx,denbx,  intlx,  int2x, 

47  $  hyx, dhyx,k1dhyx,dk1hyx,dk2hyx,k2dhyx, gamma, dgamdq, i , 
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48  $  koa123,rtsijnx,zero,q1,q2,sumx,surfno,dqi j,dqi jdz,sqng, 

49  S  dnLinbA,dhux,dhlx,e13x,cneg,cldqzl,cldqzni,cigama,koawav,tthd, 

50  +  tacoef .dacoef 

51 

52  parameterCdowni =1 .d-3,downr=1 .d- 3/0. 43429448 1 90325 18d0, 

53  +  pi=3.141592653589793238462643d0, 

54  ♦  i=(0.0d0,1.0d0),tthd=(2.d0/3.d0)»t, 

55  ♦  cneg=(0.0d0,3.141592653589793238462643d0),e13x=cneg/3.d0) 

56 

57  c  n»xlayr=maximLiii  number  of  layers  allowed 

58  c  tnxmode=maximiJii  nimber  of  modes  allowed 

59 

60  c 

61  c  use  include  file  for  parameters  of 

62  c  use  include  file  for  parameters  of 

63  c  mxlayr  max  #  layers 

64  c  mxmode  max  U  modes 

65  c 

66  Sinclude:  ‘mlaparm. i nc ' 

•••••  Begin  listing  of:  mlaparm. inc 

1  c 

2  c  include  file  to  define  the 

3  c  maximum  #  of  layers  (mxlayr) 

4  c  maximum  U  of  modes  (mxmode) 

5  c 

6  parameter  (nvxlayr=35  ) 

7  parameter  (mxmode=127) 

•••••  gnd  listing  of:  mlaparm. inc 

67  c 

68  c 

69  c****» 

70  c  acoefx-two  dimensional  complex  array  used  for  combining  two 

71  c  independent  solutions  to  stokes'  equation 

72  c  bcoefx-two  dimensional  complex  array  used  for  normalizing  height 

73  c  gain 

74  c  cqij-two  dimensional  array  containing  coefficients  for  evaluating 

75  c  qi j  in  terms  of  q11 

76  c  dqij-array  containing  coefficients  for  evaluating  qij  in  terms  of 

77  c  q11 

78  c  dqijdz-array  containing  derivatives  of  qi(z)  in  the  different 

79  c  layers 

80  c  zi-array  containing  input  hesights  for  the  modified  refractivity 

81 

82  dimension  acoefx(mxlayr, mxmode), 

83  $  bcoefxfmxlayr, mxmode), 

84  S  dqi  j{nixlayr),cqi  j(mxlayr,2),dqi  jdz(mxlayr),zi(mxlayr+1) 

85  c»***» 
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S6 

87  common  /com1/freq,waveno,sqng 

88  common  /com2/cqi j ,dqi j ,dqi jdz.nzlayr 

89  common  /papl/nrmode, koa123, surfno, zi 

90  common  /pap2/acoefx,bcoefx 

91 

92 

93  c  check  for  single  layer 

94  c 

95  c  set  a  complex  variable  koanav=- i*koa123/(waveno*waveno)  to 

96  c  avoid  repeating  computations 

97 

98  koawav=- i*koa123/(waveno*Haveno) 

99 

100  ifCnzlayr  .eq.  Dthen 

101  q1=cqi j { 1 , 1 )+zero*dqi j{ 1 ) 

102  call  surfCql , gamma, dgamdq) 

103  call  xcdai ( -q1 ,k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

104  dh2xq1 =dh2xq1+e13x 

105  int1x=cdlog{koauav*dgatndq-q1/dqi  jdz(1  ))<-2.0d0*h2xq1 

106  int2xi2.0d0*dh2xq1 -cdlogC -dqi jdz{1)) 

107  call  xcaddlsumx, intlx, int2x) 

108  rtstinx=0.5d0*sumx 

109  bcoef x( 1 ,m)=- rtsumx 

110  return 

1 1 1  end  i f 

112 

113  cldqzl=cdlog(-dqi jdz(1 )) 

114 

115  c  if  I  equals  one  then  initialize  cumulants  and  caculate  a's  and 

116  c  b's  in  bottom  layer  using  ground  boundary  conditions. 

117 

118  q1=cqij(1,1)+zero*dqi j(1) 

119  call  xcdai ( -q1 ,k2xq1 ,dk2xq1 ,k1xq1 ,dk1xql ,h2xq1 ,dh2xq1 ) 

120  dk2xq1=dk2xq1+cneg 

121  dk1xq1=dk1xq1-e13x 

122  call  surf (q1 .gamma, dgamdq) 

123  cigama=cdlog(  i»garn7ia) 

124  call  xcaddCnixnax.cldqzl -cneg+dk2xq1,cigama+cneg+k2xq1 ) 

125  call  xcadd(denax,cigama+k1xq1,cldqzl+dk1xq1> 

126  acoefxl  1  ,m)=nijnax-denax 

127  call  xcaddCdenbx, k2xq1 , acoefxd ,m)+k1xq1 ) 

128  bcoef x( 1 ,m)=-denbx 

129 

130  c  calculate  contributions  to  normalizing  integrals. 

131 

132  call  xcaddChyx, k2xq1 , acoef x( 1 ,m)+k1xq1 ) 
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133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166  c 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176  9010 

177 

178 

179  c 


hyx=bcoef x( 1 ,m)+hyx 

call  xcadd(dhyx,dl(2xq1  .acoefxd ,in)+dk1xq1 ) 
dhyx=bcoef  x(  1  ,iii)-<-dhyx 

int1x=cdlog(koawav*dgatndq-q1/dqi  jdzd  ))+2.0d0*hyx 
i nt2x=2 . 0d0*dhyx- c I dqz I 
call  xcaddlsunx, intlx, int2x) 

do  9010  l=2,nzlayr-1 
lm1=l-1 

cldqzl=cdlog(-dqi jdz(l)) 
cldqztn=cdlog(dqi  jdz(  Iml )) 
q1=cqi j( I , 1 )+zero*dqi j( 1) 

call  xcdai {-q1 ,k2xq1 ,dk2xq1 ,k1xq1 ,dk1xql ,h2xq1 ,dh2xq1 ) 

dk2xq1 =dk2xq1 +cneg 

dk1xq1=dk1xq1-e13x 

q2=cqi j ( Iml ,2)+zero*dqi j ( Iml ) 

call  xcdai ( -q2, k2xq2,dk2xq2, k1xq2,dk1xq2,h2xq2,dh2xq2) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2-el3x 

call  xcadd(hyx,k2xq2, acoef x( Iml ,m)+k1xq2) 

cal  I  xcadd(dhyx,dk2xq2,acoefx( Iml ,m)+dk1xq2) 

k1dhyx=k1xq1+dhyx 

dk1hyx=dk1xq1+hyx 

dk2hyx=dk2xq1+hyx 

k2dhyx=k2xq1+dhyx 

call  xcadd(denax,cldqznnk1dhyx,cldqzl*dk1hyx) 

call  xcadd(nunax,cldqzl-cneg+dk2hyx,cldqzm*cne9»k2dhyx) 

acoefxC  I  ,m)=rnjmax-denax 

call  xcaddldenbx, k2xq1 , acoefxl I ,m)»k1xq1 ) 

nijribx=bcoef  x(  Iml  ,m)+hyx 

dnuiibx=bcoefx(  Iml  ,m)+dhyx 

bcoefxC I ,m)=numbx-denbx 

calculate  contribution  to  normalizing  integrals. 

int 1x=cdlog( -ql/dqi jdz( I )+q2/dqi jdzl Iml ) )+2 . 0d0*numbx 
call  xcaddCsumx.sumx, intlx) 
call  xcadd(dhux,dk2xq1 , acoef x( I ,m)+dk1 xql ) 
dhux=bcoefx( I ,m)+dhux 
i  nt  1  x=2 . 0d0*dnLitfcx-c  Idqzm 
i nt 2x=2 . 0d0*dhux - c I dqz I 
call  xcaddlsumx.sunx, intlx) 
r.  xcaddCsmw.sunx,  int2x) 
continue 


if  I  equals  nzlayer,  calculate  a's  and  b's  using  outgoing 
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184 

185 

186 

187 
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189 

190 
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195 
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200 
201 
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204 
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207 

208 

209 

210 
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213 

214 
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216 

217 

218 

219 

220 
221 
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223 

224 

225 

226 


c  wave  in  top  layer. 

nzm1=nzlayr- 1 

q1=cqi j (nzlayr, 1 )+zero*dqi j(nzlayr) 

call  xcdai ( -q1 , k2xq1 ,dk2xq1 , klxql ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

dh2xq1=dh2xq1+e13x 

q2=cqi j  Cnzmi , 2)+zero*dqi j (nzral ) 

call  xcda i ( - q2 , k2xq2 , dk2xq2 , k 1 xq2 , dk 1 xq2 , h2xq2 , dh2xq2 ) 

dk2xq2=dk2xq2+cneg 

dk 1 xq2=dk 1 xq2 - e 1 3x 

cal  I  xcadd(hyx,k2xq2,acoefx(nzm1  ,m)'rk1xq2) 
niiit>x=bcoefx(nzlayr- 1  ,m)+hyx 
bcoefx(nzlayr,iii)=nLfTibx-h2xq1 

c  calculate  contribution  to  cumulants. 

int 1x=cdlog( -q1 /dqi jdzCnzlayr )*q2/dqi jdzlnzml ) )♦ 

S  2.0d0*numbx 

call  xcadd(sutnx,s'jmx,  intlx) 

cal  I  xcadd(dhyx,dk2xq2,acoefx(nzm1 ,m)+dk1xq2) 
dnumbx^bcoef  x(nziTil  ,m)+dhyx 
int  1x=2.0d0*dnu(tibx-cdlog{dqi  jdz(nzm1 ) ) 
call  xcadd<sunx,sumx,  intlx) 
dhux=bcoefx(nzlayr,m)+dh2xq1 
int2x=2.0d0*dhux-cdlog( -dqi jdzCnzlayr) ) 
call  xcaddCSLTix.sumx,  int2x) 

c  renormalize  b's  so  that  height  gain  integral  equals  un’ty. 

rtsijnx=  .5d0*sumx 
do  9900  1 1  =  1 , nzlayr 

bcoef x( 1 1 ,m)=bcoef x( 1 1 ,m) - rtsumx 
9000  continue 

^***************************»*********************************** 


l=nzlayr 
Im1 =1-1 

cldqzm=cdlog(dqi jdz( Iml )) 
c Idqz l=cdlog( -dqi jdz( I ) ) 

c  calculate  q  and  associated  quantities  at  bottom  of  layer  I 
q1=cqi jCl,1)+zero*dqi j(l) 

call  xcdai ( -q1 , k2xq1 ,dk2xq1 , klxql ,dk1xq1 ,h2xq1 ,dh2xq1 ) 
dh2xq1=dh2xq1+e13x 
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227 

228 

229 

230 
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232 

233 
23A 
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236 

237 

238 

239 
2^10 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 
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253 

254 

255 

256 

257 

258 

259 

260 
261 
262 

263 

264 

265 

266 

267 

268 

269 

270 

271 

272 

273 


q2=cqi  j  ( Itnl  ,2)+zero*dqi  j  ( Iml ) 

ca 1 1  xcda i  ( ■  q2 , k2xq2 , dk2xq2 , k1xq2, dkl xq2 , h2xq2 , dh2xq2 ) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2-e13x 


c  Caculate  acoefxC  Iml  ,m)  ,bcoefx(  litil  ,m) 

c  and  cunulants  using  outgoing  wave  in  nzlayr 

^***** 

dh2k1x=dh2xq1+k1xq2 

h2dk1x=h2xq1+dk1xq2 

h2dk2x=h2xq1+dk2xq2 

dh2k2x=dh2xq1+k2xq2 

cal  I  xcadd(denax,cldqzt-cneg+dh2k1x,cldqzm+cneg+h2dk1x) 
cal  I  xcadd(numax,cldqzm+h2dk2x,cldqzl+dh2k2x) 

c  If  in  the  nzlayr-1  layer  the  magnitudes  of  A  coefficients  from 
c  integration  up  and  down  differ  by  less  than  0.02  dB  and  their 
c  phases  differ  by  less  than  O.OOlpi,  the  A  and  B  coefficients 
c  obtained  from  integration  up  will  be  accepted. 

tacoef =numax-denax 
dacoef=tacoef -acoefxC Iml ,m) 
dif r=dabs(dreal (dacoef )) 
if  (difr  .It.  downr)  then 
dif  i=ditnag(dacoef  )/pi 
di f i =dabs<di f i -dni nt (di f i /2.d0)’2.d0) 
if  (difi  .It.  downi)  return 
end  i  f 

acoefxC Iml ,m)=tacoef 

call  xcaddCdenbx, k2xq2, acoefxC Iml ,m)+k1xq2) 
bcoefxC Iml ,m)=h2xq1 -denbx 

c  calculate  contributions  to  cumulants 

SLiTw=cdlogC  -ql/dqi  jdzc  I  )+q2/dqi  jdzC  Iml)  )+2.0d0*h2xq1 

call  xcaddCdh I x,dk2xq2, acoefxC Iml ,m)*dk1xqZ ) 

dhlx=bcoefxC Iml ,m)*dhlx 

int1x=2.0d0*dh2xq1 -cldqzl 

call  xcaddC intlx.sumx, intlx) 

i nt  2x=2 . 0d0*dh I x - c I dqzm 

call  xcaddCsimx, intlx, int2x) 

do  9030  l=nzlayr-1 ,2,-1 
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27A 

275 

276 

277 

278 

279 

280 
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287 
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290 
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308 
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311 

312 

313 

314 

315 

316 

317 

318 

319 

320 


lm1=l-1 

cldq2l=cdlog(-dqi  jdzCD) 
cldqzm=cdlog(dqi jdz( Iml ) ) 

c  calculate  q  and  associated  quantities  at  bottom  of  layer  I 
q1=cqi j( I , 1 )+zero*dqi j ( I ) 

call  xcdai(-q1 ,k2xq1 ,dk2xq1,k1xq1 ,dk1xq1 ,h2xql,dh2xq1 ) 
dk2xq1 =dk2xq1+cneg 
dk1xq1=dk1xq1 -e13x 

q2=cqi j( Iml ,2)+zero*dqi j( Iml ) 

cal  I  xedai (-q2,k2xq2,dk2xq2,k1xq2,dk1xq2,h2xq2,dh2xq2) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2-e13x 

dh2xq2=dh2xq2'»el3x 

c  Calculate  acoefx( Iml ,m),bcoefx( Iml ,m)  and  cumulants 

c  using  continuity  relations  in  terms  of  the  linearly 

c  independent  functions  k1  and  k2 

call  xcadd(hyx, k2xq1 ,acoefx( I ,m)*k1xq1 ) 

cal  I  xcadd(dhyx,dk2xq1 ,acoefx( I ,m)tdk1xq1 ) 

k1dhyx=k1xq2+dhyx 

dk1hyx=dk1xq2+hyx 

dk2hyx=dk2xq2+hyx 

k2dhyx=k2xq2+dhyx 

call  xcaddldenax.cldqzl ■cneg*k1dhyx,cldq2tTH-cneg»dk1hyx) 

cal  I  xcadd(numax,cldqzm*dk2hyx,cldqzl*k2dhyx) 

acoefxC  Iml  ,m)=nun',ax-denax 

call  xcaddldenbx, k2xq2,acoofx( Iml ,m)+k1xq2) 

numbx=bcoef x( I ,m)*hyx 

dnijiibx=bcoef  x<  I  ,m)*dhyx 

bcoefxC  Iml  ,m)=numbx-dcn.bx 

c  calculate  contributions  to  cumulants. 

int  1x=:cdlog(  -ql/dqi  jdz(  I  )+q2/dqi  jdz(  Iml )  )  +  2.0dO*numbx 

call  xcaddCsLrnx.sumx,  intlx) 

call  xcadd(dhlx,dk2xq2,acoefx( Im1,m)+dk1xq2) 

dhlx=bcoefx( Iml ,m)+dhlx 

i nt 1 x=2 . OdO*dnumbx - c I dqz I 

int2x=2.0d0*dhlx-cldqzm 

call  xcadd(SLnx,sunx,  intlx) 

call  xcaddlsumx, sumx, int2x) 
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9030 


cont i nue 


321 

322 

323 

324 

325 

326 

327 

328 

329 

330 

331 

332 

333 
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335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 


c  if  I  equal  to  one  calculate  ground 
c  contribution  to  cuntulants  and  renormalize  bcoefx's 

1  =  1 

q1=cqi j ( I , 1 )+zero*dqi j ( I ) 

cal  I  xcdai ( -q1 , k2xq1 ,dk2xq1 , klxql ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

dk2xq1=dk2xq1+cneg 

dk1xq1=dk1xq1 -e13x 

call  xcaddChyx, k2xq1 , acoefxC I ,m)»k1xq1 ) 
cal  I  xcadd(dhyx,dk2xq1 ,acoefx( I .ml+dkixql ) 
call  surf (q1 , gamma, dgamdq) 
numbx=bcoefx( I ,m)+hyx 
dnLribx=bcoef  x(  I  ,m)+dhyx 

int1x=cdlog(koawav*dgamdq-q1/dqi jdzl I ))f2.0d0*numbx 
i  nt2x=2.0d0*dn(jnbx-cdlog(  'dqi  jdz(  1 ) ) 
call  xcadd(sumx,suinx,  intlx) 
call  xcaddCsunx.sumx, int2x) 

c  renormalize  b's  so  that  height  gain  integrals  equal  unity. 

rtsumx= .5d0*sumx 

do  9020  1 1  =  1 ,nzlayr-1 

bcoefxC 1 1 ,m)=bcoefx( 1 1 ,m)-rtsumx 
9020  continue 

bcoefx(nzlayr,m)=- rtsamx 


return 

end 
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